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Executive Summary 


This report summarizes our derivations of analytical expressions for the frequencies and damping constants 
for small-amplitude axisymmetric shape oscillations of a liquid drop suspended in an immiscible fluid host in 
microgravity. In particular, this work addresses large Reynolds number shape oscillations and focuses on the 
surface rheological effects that arise from the presence of insoluble surfactants at the interface. Parameters 
characterizing viscous effects from the bulk phases, surface viscous effects, Marangoni effects from the surface 
advection and diffusion of surfactants, and the Gibbs elasticity are all considered and analyzed to determine 
the relative importance of each contribution. 

Supplementing the analytical treatment for small-amplitude oscillations, a numerical boundary integral 
equation formulation is developed for the study of large-amplitude axisymmetric oscillations of a drop in 
vacuum. The boundary integral formulation is an extension of classical potential flow theory and approx- 
imately accounts for viscous effects in the bulk fluid as well as the surface viscous and Marangoni effects 
resulting from an insoluble surfactant contaminating the interface. 

Theoretical and numerical results are presented for four distinct cases. These range from the case when the 
effects of the surfactants are “negligible” to “large” when compared to the viscous effects in the bulk phases. 
The feasibility of the non-contact measurement of the surface parameters, using experimental observations 
for the oscillation frequencies and damping constants of drops and bubbles, is discussed. 
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Chapter 1 

Introduction 


The dynamics of the shape deformations of drops and bubbles are important in a wide variety of prac- 
tical applications including heat-transfer [7], chemical engineering [16], multiphase flow [24], and nuclear 
physics [13]. Any system in which the surface area to volume ratios are high may be greatly influenced 
by the presence of surfactants. Surfactants are long molecules with separated hydrophilic and hydrophobic 
segments that preferentially adsorb to an interface and give it viscoelastic properties. For a Newtonian inter- 
facial fluid these viscoelastic properties are characterized by surface tension, a Gibbs elasticity, and surface 
shear and dilatational viscosities. Since the dynamics of drops and bubbles depend on these viscoelastic 
properties, experiments have been proposed for the non-contact measurement of these properties through 
the observations of forced and freely oscillating drops and bubbles [48, 20, 9, 25, 15]. The present work aims 
to determine analytical expressions for the dependence of the frequency and damping constants on these 
viscoelastic interfacial properties, needed for extracting their values from the experimental data. 

The study of the small-amplitude shape oscillations of drops and bubbles has a long history. By using 
a normal-mode analysis to separate the time variable, Rayleigh [39] determined the small-amplitude linear 
frequencies of an inviscid drop oscillating in vacuum. This analysis was generalized to treat the case of an 
inviscid drop submerged in an inviscid medium by Lamb [23], who also used the inviscid velocity solutions 
to estimate the rate of damping of a weakly viscous drop in vacuum or a gas bubble immersed in a weakly 
viscous medium. A similar dissipation estimate was obtained for arbitrary viscosities by Reid [40] for a 
drop in a vacuum or low density gas, and by Valentine et al. [53] for a drop submerged in a fluid medium. 
Miller and Scriven [29] used a normal-mode analysis of the linearized Navier-Stokes equations to identify the 
primary role of boundary layer dissipation on the frequency shift and damping rate for the weakly viscous 
drop arid medium system. Their formulation also included viscoelastic surface properties, but they limited 
their analysis to the case of free and inextensible interfaces. Prosperetti [37, 38] used a more general technique 
based on Laplace transforms to study the transient regimes of a weakly viscous drop and medium system 
and showed that the normal-mode results are recovered in the asymptotic limit for long times. Higher-order 
corrections to Miller and Scriven ’s normal-mode analysis have been obtained by Marston [28] and Asaki and 
Marston [6] for free and acoustically forced shape oscillations of a weakly viscous drop and medium system. 
Lu and Apfel [25] have extended Miller and Scriven’s normal-mode analysis to include the effects of a soluble 
surfactant in the outer medium. A similar normal-mode analysis has been used by Tian et al. [45] to analyze 
the effects of a soluble surfactant on a drop oscillating in vacuum in the quadrupole mode. 

Due to its relative complexity, the large-amplitude analysis of drops and bubbles has been limited to 
inviscid or numerical models. Tsamapoulos and Brown [51, 52] and Natarajan and Brown [33, 34] have used 
inviscid models to study nonlinear shape oscillations and mode coupling. Foote [19] and Alonzo [2] analyzed 
the nonlinear oscillations of a viscous drop in vacuum using the marker-and-cell numerical method. The large- 
amplitude oscillations of a viscous drop has also been studied by Basaran [7] using a finite-element method. 
The boundary integral numerical method for potential flow has been used by Lundgren and Mansour [26] 
and Shi and Apfel [44] to simulate high Reynolds number drop oscillations, where weak viscous effects were 
included in the formulation by a modification of the normal stress boundary condition. The weakly viscous 
analysis using the boundary integral numerical method has been extended for the case of charged drop 
forced in an electric field by Feng and Beard [18] and an acoustically forced drop incorporating the effects 
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of constant surface viscoelastic properties by Chen et al [15]. 

The shape oscillation experiments of an acoustically levitated drop or bubble began with Marst.on and 
Apfel [27]. Trinh et al [46, 47, 48] have also used acoustic forcing to study the small and large-amplitude 
shape oscillations of a neutrally buoyant drop submerged in an outer medium. Asaki and Marston [6, 5] have 
analyzed acoustically positioned gas bubbles oscillating in the presence of insoluble and soluble surfactants 
and Trinh et al [49] have studied the combined effects of acoustic and electric field forcing on charged drops 
in air. Recently Apfel et al [3] have studied large- amplitude oscillations for a surfactant contaminated drop 
in air in microgravity. 

The present work considers high Reynolds number shape oscillations of a drop and medium system and 
extends the knowledge of small-amplitude shape oscillations in the presence of surfactants with the derivation 
and analysis of a total mechanical energy equation for the system. The analysis of this total mechanical energy 
equation with matched asymptotic expansion techniques generalizes the energy approaches of Lamb [23] and 
Hsu and Apfel [20] for the determination of the oscillation frequency and damping constants for the system 
and provides physical insight into the origin and significance of resulting terms. The understanding of large- 
amplitude shape oscillations of a drop in vacuum is also improved in the present work with the development 
and analysis of an efficient boundary integral numerical method for high Reynolds number shape oscillations 
that incorporates, in an approximate way, the viscous effects in the bulk fluid and surface viscoelastic effects 
arising from the presence of an insoluble surfactant at the interface. 
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Chapter 2 

Theory 


2.1 Motivating example: An oscillating flat plate 

Figure 2.1 shows a flat plate of mass m submerged in a fluid of density p and viscosity p and attached to 
a wall by a spring of stiffness k. The equations and boundary conditions describing the one-dimensional 
time-dependent displacement of the plate x(f) from equilibrium and fluid velocity u(y,t) are given by 



Du 

rnx(t) + kx(t) = .4/x— (0, t ) , 

(2.1) 


a:(0) = x 0 , i(0) = 0 , 

(2.2) 

and 

Du D 2 u 

(2.3) 


u(y , 0) = 0 , u( 0, t) = x(t ) , u(oo, t) — 0 . 

(2.4) 

Gravity and end effects have been neglected. Equation (2.2) corresponds to the plate being released from 
an initial displacement x Q . The right-hand side of (2.1) represents the total viscous force on the plate with 
total surface area A. The above equations and boundary conditions may be nondimensionalized with the 
time scale lj" 1 = ( k/m ) 1 / 2 and the length scale x G . The distance normal to the plate is scaled with the 
Stokes’ boundary layer thickness (/// puJ Q ) 1 and the nondimensional form of the equations and boundary 


conditions are 



Figure 2.1: a) Schematic of oscillating plate system, b) Stokes 7 boundary layers above plate. 
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(2.5) 


and 


X{t) +x(t) — £ — (0J.) , 

oy 

x(0) = 1 , ±(0) = 0, (2.6) 


du 

dt 


d 2 u , 
dxj 2 


(».0 = iri (».*). 


(2.7) 


u(^0) = 0, u(0, t) = x(t ) , u(oo,t) = 0. 


( 2 . 8 ) 


Equations (2. 5-2. 8) contain the single parameter 


e — 



(2.9) 


which may be interpreted as the ratio of the mass of the fluid in the boundary layer to the mass of the plate. 

In the following sections the motion of the plate x(t) when e < 1 is analyzed using two methods. 
Comparisons are made with shape oscillations of the drop/medium system. 


2.1.1 Exact solution using Laplace transforms 

Equations (2. 5-2.8) may be solved exactly using Laplace transforms in the time variable. The Laplace 
transform is defined such that 

x(s) = C{x(t)} = f x{t)e~ $l dt . (2.10) 

Jo 

Equations (2. 5-2. 8) and (2.10) give the solutions for the plate and fluid motion in transform space as 

tS + £\/s 


x(s;e) = n 

(1 *+* S 2 ) + ESyfs 
u{yTS\e) = (sx-l)e-^y. 


( 2 . 11 ) 

( 2 . 12 ) 


When e = 0, the motions of the plate and fluid are time harmonic with the fluid motion above the plate 
corresponding to Stokes’ 2nd Problem: 

x(t) = £ _1 {:r(s;0)} = cos (t) (2.13) 

u{y,t) = £ _l {LZ( 2 /, s\ 0)} = e~ y ^ cos(t — y/V2 ) . (2-14) 


If ^ > 0, the time-dependent motion of the plate may be expressed in terms of an inverse Laplace transform 


x(t) = C 1 {a:(«;e)} 


(s + £^fs)e st ds 

2 7T2 J c _ loc (1 + S 2 ) +ESy/s 


(2.15) 


and evaluated in the complex 5 -plane using residue theory. See Figure 2.2 for a description of the integration 
contour. Note that for small £, the poles of the integrand in (2.15) are shifted from ±i by terms of O(e) 
which are readily calculated. With a total solution of the form x(t) = x sp (t.) + £b r (0> the leading-order 
contribution from the two simple poles is 


x sp (t) = e- £t / 2 ^cos[(l - e/2y/2)t] + O(ser^) . 


(2.16) 


Here x sp (t) is an exponentially damped time harmonic oscillation with a frequency slightly shifted from 
unity. The leading-oder contribution from the branch cut 


, x -e f°° Vte~ rt 

Xhr ® ~ n Jo (1 + r 2 ) 2 + e 2 r 3 T 

(2.17) 

-el , n 

~ — == as t — ► oo ; e — > (J . 

2y/n t 

(2.18) 
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Figure 2.2: Contours around two simple poles and a branch cut in the complex s-plane. 

Here x^ r (t) provides a small correction that decays algebraically for long times. Its contribution dominates 
that from the simple poles for long times, but by that time the motion of the plate is negligibly small. 

Similar behavior has been seen in the analysis of the shape oscillations of drops and bubbles using Laplace 
transforms. Prosperetti [37, 38] has shown that the time dependent amplitude of the shape oscillations of 
drop/medium system is composed of a discrete spectrum of exponentially decaying functions and a small 
continuous spectrum. Roberts and Wu [41] have shown that the continuous spectrum for the shape oscilla- 
tions of a bubble has an algebraic decay and therefore dominates the total solution for long times. Each of 
these treatments contains a level of analysis that is impractical when applied to the much more complicated 
case of the shape oscillations of drops and bubbles when surfactants are present. The following section 
outlines an averaging method which allows for the calculation of the leading order solutions in a simple way. 


2.1.2 Approximate solution using an averaging method 

The analysis using the averaging method [21] begins with the derivation of a total mechanical energy equation 
for the system. For the plate and fluid system this involves multiplying the plate equation (2.5) by x{t) and 
the fluid equation (2.7) by u(y,t) and integrating over y to yield 


d 

dt 


r x 2 x 2 " 1 

~2 + T 


= £X^(0,t) 

dy 


(2.19) 


and 


r °° u 2 


dt 


L \ dy = 


dru 
1 dy 2 


^ d y 


= -±^(0,t) 

dy 


f(£) v (220) 

These two equations may now be combined to form the total mechanical energy equation for the system 


d 

dt 




(2.21) 


The left-hand side of (2.21) is the time rate of change of the total energy of the system, which consists of the 
kinetic energy of the plate and fluid and the potential energy in the spring. The right-hand side of (2.21) 
is the the total viscous dissipation in the fluid. For e = 0, this equation shows that, the total energy in the 
system is conserved. 

The averaging method uses the time-periodic solutions of the plate and fluid to find the time-average 
of the quantities of kinetic energy, potential energy, and dissipation rate. An energy equation is then 
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constructed which yields terms with equivalent time averages. This new energy equation is then simplified 
to form a linear oscillator equation that captures the leading order behavior of the plate. This procedure is 
demonstrated below. 

Suppose in general that x(t) = SijC/e 1 *} and u(y,t) = $1{U exp[— (1 + i)y/y/2 + i£]}, where U is a complex 
amplitude. (For the particular boundary conditions above U == 1.) It is easily shown that 



( 2 . 22 ) 

(2.23) 


Here the angular brackets {) denote the time average over one period of oscillation. Replacing the fluid 
terms in (2.21) with terms involving x on the right-hand sides of (2.22) and (2.23), which have the same 
time average, yields 


d_ 

dt 


+ £ 


x 2 x 2 

2V2 + 2 


= -s 


v/2' 


(2.24) 


Taking the time derivative of the left-hand side of (2.24) and dividing through by x(t) gives the following 
oscillator equation 

(! + + 7f i: + a: = 0, t 2 ' 25 ) 

which contains added mass and damping terms. The leading order solution is 


x(t) = e _et/2 ' /5 cos[(l - E/2y/2)t] + 0{e 2 ) . 


(2.26) 


Comparing equations (2.16) to (2.26) shows that the averaging method captures the leading order contribu- 
tions to the time dependent motion of the plate, but fails to capture the additional O(e) corrections from 
the branch cut. 


2.2 Drop oscillations: equations and boundary conditions 

Returning to the problem of drop oscillations, consider a liquid drop of density p and viscosity p suspended 
in an infinite fluid medium of density p and viscosity //, in the absence of gravity. Both fluids are assumed to 
be incompressible and Newtonian. The continuity and momentum equations in each phase take the forms 


<1 

< 

II 

o 

p—~ — V ■ 
h Dt 

n 

for x € V m (t) , 

(2.27) 

<1 

II 

o 

t> 

ii 

n 

for x £ V m (t ) . 

(2.28) 


Here v and v, respectively, refer to the drop and medium velocity fields, V m (t) and V m (£) are the material 
volumes of the drop and medium, and the stress tensors II and II are given by 

n=-pI + 2/iE, E = l[(Vv) + (Vv) T ] , (2.29) 

n = —pi + 2fiE , E = h(Vv) + (Vv) T ]. (2.30) 

In (2.29) and (2.30), p and p represent the pressures in the two fluids, I is the isotropic unit tensor, and E 
and E are the symmetric and traceless rate-of-strain tensors. 

These field equations need to be supplemented by boundary conditions at infinity and at the material 
interface S m (t ) between the drop and medium. At infinity, the velocity field vanishes and the pressure tends 
to a constant value. The interface is assumed to be covered with surfactants and therefore possesses its own 
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rheological properties, which may be characterized by the surface stress tensor II S [16]. The no-slip and 
stress balance boundary conditions at the interface thus assume the respective forms [16] 


v 

n (n - n) 


* i 

-v s n s / 


for x € S m (t) . 


(2.31) 


The velocity of the interface is equal to the fluid velocity in the drop or medium evaluated at S m (t ). For 
convenience, this velocity is subsequently denoted by v. The surface stress tensor is also assumed to be 
“Newtonian” and defined by a Boussinesq-Scriven constitutive relationship of the form [16, 31, 43] 

n s = al s + + K,I S (V S • V) . (2.32) 


Here cr, /x s , and k s , respectively, refer to interfacial tension, surface shear viscosity, and surface dilatational 
viscosity. Also, I s = I - nn is the surface unit tensor and n is a unit vector normal to the interface pointing 
into the surrounding fluid medium. V s = I s • V is the surface gradient, and E s is the symmetric and traceless 
surface rate-of-strain tensor defined by 

E s = l[(V.v) I s + I s (V s v) 7 '] - |l.(V, • v) . (2.33) 

In this work (j s and k s are taken to be constant. The surface tension cr, however, depends on local 
surfactant concentration T. For small-amplitude drop oscillations the concentration of surfactants T is 
assumed to vary only slightly from the equilibrium concentration F eq and the surface tension is approximated 
as a linearly decreasing function of the surfactant concentration [16] 

a(r)=a eq -^(r-T eq ). (2.34) 

1 eq 

Here a eq is a constant equilibrium surface tension and e s is the Gibbs elasticity. 

The surfactant transport equation for an insoluble surfactant is given by [31] 

^ + V s • (vr) = D S V S 2 r , (2.35) 

at 

where D s is the surface diffusivity of surfactants and v is the velocity of the interface. The surfactant 
transport equation is coupled to the equations governing the bulk fluid motions through (2.34). 


2.3 The total mechanical energy equation 


The total mechanical energy equation is obtained by dot multiplying the momentum equation in (2.27) 
by v, the momentum equation in (2.28) by v, integrating over the respective material volumes V m (t) and 
and adding the resulting equations. With the aid of the bulk Reynolds Transport and Divergence 
Theorems [4, 31], the no-slip boundary condition (2.31), and the incompressibility conditions (2.27) and 
(2.28), the following equation is obtained 


d 

dt 



-p\v\ 2 dV + 


L> fdv }= 

- 2p(E : E) dV - J 2/t(E : E) dV 
-jT h-(tl-U)-vdS. 


(2.36) 


The integrand in the last term on the right-hand side of (2.36) may be simplified using the stress balance 
boundary condition (2.31): 


n (ri-n) v 


(-v s -n s )-v 

nj : V s v - V s • (IC • v) 

2p s (E s : E s ) + k s (V s -v) 2 
+rr(V s • v) - V, • (n s • v) , 


(2.37) 
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so that 


c i_ 
dt 


\p\v\ 1 dV + 

V'm 


IJ mUv } = 


2/i(E : E) dV 


i\„, 




2ft(E : E) dV 


-f 2fi s (E s : E s ) dS - f k s (V s -v) 2 

J S n i J Sm 

-j <t(V s • v) dS - j s V s • (n s -\)dS. 


dS 


(2.38) 


The last term on the right-hand side of this expression is zero by the the Surface Divergence Theorem [31], and 
the second to last, term may be simplified using (2.34) and the Surface Reynolds Transport Theorem [31]. 
The general forms of the Surface Reynolds Transport and Surface Divergence Theorems are given in the 
Appendix A. The final form for the total mechanical energy equation is 



(2.39) 


The terms on the left-hand side represent the time rate of change of the total kinetic energies in each phase 
and the potential energy. The first two terms on the right-hand side of (2.39) represent the viscous dissipation 
rate in the bulk phases. The next two terms are similarly identified as dissipation terms arising from surface 
shear and dilatational viscosities. The last term on the right-hand side of (2.39) contains the surface tension 
gradient, or Marangoni, effects and couples the total mechanical energy equation to the surfactant transport 
equation. Its interpretation is less obvious and is discussed below. 

If the bulk and surface viscosities and the Gibbs elasticity e s are set to zero, then equation 

(2.39) shows that the total energy of the system is conserved. If only the Gibbs elasticity is set to zero, total 
viscous dissipation rate has contributions from each of the bulk phases with additional contributions from 
the surface phase. 

An interpretation of the Marangoni term may be seen through an analysis of the surfactant transport, 
equation (2.35), which may be equivalently expressed in the following form 


r\ 

r e ,(v,-v) = --(r-r e ,)-v v s (r-r e? ) 

-(r - r e ,)(v, • v) + d s v s 2 (r - r„) . 


(2.40) 


After multiplying this equation by e s (r - Y eq )/T' 2 e<r integrating over the material surface S ln {t), and using 
the Surface Reynolds Transport and Surface Divergence Theorems, the Marangoni term in (2.39) may be 
re-expressed as 

^,r-r„ K v,.v„5 - -±j^-T.,)Us 

- Is ^|v s (r-r e ,)| 2 ds 

-f r^(r-r cg ) 2 (v., -v) dS. (2.4i) 

The Marangoni term must therefore be interpreted as having several contributions. The first contribution is 
the time rate of change of an additional energy storage term arising from the non-equilibrium distribution of 
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surfactants. The second contribution is dissipative and arises from the irreversible diffusion of surfactants. 
The last term on the right-hand side of (2.41) is not clearly identifiable as a dissipative or energy storage 
term. For small-amplitude oscillations this term will be shown to be negligibly small compared to the other 
two. Depending on the surface Peclet number characterizing the surfactant transport, the leading-order 
behavior of the Marangoni term can be shown to be dissipative, provide an additional energy storage, or 
give a combination of these [32]. 


2.4 Analysis of total mechanical energy equation 

In this section the averaging method is used to analyze the total mechanical energy equation of the drop/medium 
system performing small-amplitude shape oscillations. 

Nondimensionalization is based on the inertial scales from the drop fluid properties: 

/ pH 3 \ l / 2 

length = R, time = cj -1 = ^ j , mass — pR 3 . (2.42) 

Surface tension and surfactant concentration are scaled on their equilibrium values a eq and T e(r The nondi- 
mensional forms for the total mechanical energy equation and surfactant transport equation are thus 


and 


where 


d 

dt 

{K.E. 

+ P.E.} 

= 

— a 2 {Bulk Diss.} 






— /r* {Surf. Diss.} 






H-e* {Marangoni} 

(2.43) 


ar 

dt 

+ V s ■ (\ 

r F) : 

= ^v s 2 r, 
Pe s 

(2.44) 

K.E. 



dV 

+ - / hv| 2 dV 

(2.45) 



Jv m 2 1 


P ' v;„ 2 

P.E. 

= 

L ds 



(2.46) 

Bulk Diss. 



[ 2(E: 

E), 

dV + 1 - 1 2(E : E) dV 

(2.47) 



Jv m 

P Jv„ 


Surf. Diss. 

— 

1 2(E S 

■ E. 

i )dS+'^ [ (V s • v) 2 dS 

(2.48) 



Js m 


P*s Js,„ 


Marangoni 

= 

fsJ T ~ 

1)(V S -\)dS . 

(2.49) 


Equations (2.45) (2.49) represent, respectively, the total nondimensional kinetic and potential energy, the 
bulk and surface viscous dissipation rate, and the Marangoni term. 

The nondimensional Marangoni term can be re-expressed as 

{Marangoni} = — — {Mar. E.} 

{Mar. Diss.} 

a e$ 

— { Remainder } (2.50) 


where 


Mar. E. 



^(r - \) 2 ds 


(2.51) 
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Remainder 


(2.53) 


= £ (r- i)' 2 (v s 


v) dS . 


Mar. E. represents the nondimensional Marangoni stored energy term. Mar. Diss. is a dissipation rate. The 
remainder term will be shown to be negligible for small oscillations. 

In addition to /x / /i and p/p, the dimensionless parameters in (2.43) (2.53) are defined to be 


(Q 2 ,/i;,«’,e*,Pe s ) = 


/i ft S' J -' ^ s 


j P R y 


' eq 


' eq 


uiR 2 

~D~ S 


(2.54) 


These represent, respectively, the inverse of the Reynolds number, the dimensionless surface shear and 
dilatational viscosities, the Gibbs elasticity, and the surface Peclet number. 

In order to analyze the drop/medium system with the averaging method, the time-periodic velocity pro- 
files in each phase must be determined. The shape oscillations are assumed to be small and axisymmetric, 
and the Reynolds number is assumed to be large so that deviations from potential flow are confined to 
thin Stokes boundary layers near the interface. The velocity and pressure fields in each phase are there- 
fore found by the solution of a singular perturbation problem, which may be approximated using matched 
asymptotic expansion techniques. For this purpose, the nondimensional shape of the interface and surfactant 
concentration are taken to be 


r(6,t\e) - 1 + eat(t)Pt(co®0) + 0(e*) (2.55) 

T (6,t;e) = 1 +egt{t)Pi{cos6) + 0{e 2 ) . (2.56) 


Here £ is a small dimensionless parameter used to linearize the equations and Pt is the Legendre polynomial 
of order i , where 0 represents the polar angle measured from the axis of symmetry. The time-dependent 
amplitudes for the axisymmetric perturbations in shape Qf(t) and surfactant concentration gt(t) are given 
by the real parts of the complex quantities 

a e {t) = 3t{Ae’ n ' 0 *} (2.57) 

gt{t) = . (2.58) 

Here .4 and G are complex amplitudes. Since the Reynolds number is assumed to be large, the assumed 
nondimensional frequency is taken to be the base frequency for small-amplitude inviscid oscillations of 
a drop/medium system given by the Lamb formula [23] 


2 = l(l-l)(< +!)(* + 2) 

" [(*+1 )+* p / p ] 


(2.59) 


Using the approximate forms for the velocity profiles found from the asymptotic analysis, order-of- 
magnitude estimates are used to identify the dominant contributions to each term in the total mechanical 
energy equation. Time- averaging these dominant contributions over one oscillation period further simplifies 
the terms and allows for the derivation of a damped harmonic oscillator equation. The characteristic fre- 
quencies for this equation capture the leading-order behavior of the drop/medium system influenced by the 
effects of viscosity. When surfactants are introduced, the damped harmonic oscillator equation is coupled to 
the surfactant transport equation through the Marangoni term. For that case, the leading-order behavior of 
the system is described by the simultaneous solution of two coupled equations. 

The following sections begin with an analysis of the base potential flow for the shape oscillations of 
the drop/medium system. The succeeding four cases include the effects of viscosity in the system with 
increasingly greater effects from the presence of an insoluble surfactant at the interface. 


2.4.1 Base flow: inviscid oscillations 

This section describes the analysis of the inviscid drop/medium system with no surfactants. All the dissipa- 
tion terms on the right-hand side of the total mechanical energy equation (2.43) are zero and the total energy 



of the system is conserved. This case serves as the base flow for large Reynolds number shape oscillations 
and demonstrates the application of the averaging method to the drop/medium system. 

For the shape perturbation in (2.55), the inviscid flow fields in each phase are given by the solution of a 
regular perturbation problem. If the flow is irrotational, the scalar velocity potentials (scaled with ojR 2 ) for 
the flow in each phase may be expanded in regular perturbation series in the small parameter e 

(f>(r,8,t]£,a) ~ £(po(r,6,t) + 0(e 2 ) (2.60) 

(j){r,0,t\£,a) ~ £0o(r,0,t) + 0(e 2 ) , (2.61) 

Here the arbitrary constant in the scalar velocity potentials has been set to zero in both phases. These 
potentials are found by solving Laplace’s equation within each phase 

V 2 <f) o=0 for r < 1 
V 2 0o = 0 for r > 1 

The potentials are subject to the 0{£) kinematic boundary conditions 

?p. = dp = *n»Ae <n<ot P/(cos0) at r = 1 (2.64) 

or or 

and the conditions that the potential in each phase remain finite. 

The nondimensional pressure (scaled with pu) 2 R 2 ) is related to the potential in each phase through the 
linearized unsteady Bernoulli equation 

P = 

P = 

Here p eq and p eq are the equilibrium pressures in each phase satisfying the nondimensional \oung-Laplace 
equation p eq - p eq (p/p) = 2. 


Peq -e^ + 0(e 2 ) (2.65) 

~ P “ £ ^w) + 0{£2) ■ {2M) 


(2.62) 

(2.63) 


Time-periodic solutions 

For at{t) = 3£{Ae* n£0 *}, the potential flow equations may be readily solved. The resulting velocity compo- 


nents v = V 0 — e r v r 4 - e$ve and pressure in the drop phase are the real parts of 

v r o(r, 0 ,t;e) = £iQ(oAe’ ntot Pe(cos 0 )r e ~ 1 + 0 (e' 2 ) ( 2 . 67 ) 

v M = etn to Ae <n «‘l^(cosfl)r , - l + 0 (e 2 ) (2.68) 

p(r, 9 ,t-,e) — Peq - en^Ae^jPficosO)^ + 0 {e 2 ) ( 2 . 69 ) 

The velocity components v — V <5 = e r v r + gqVq arid pressure in the medium phase are similarly 

t > r0 (r,0,t;e) = ein (0 Ae mtot P ( {cos6)r- {l+2) + (D(e 2 ) ( 2 . 70 ) 

veo (r,9,t;e) = -«n w Ae m '»‘^ T y^(cosfl)r-« /+2 > + 0(e 2 ) ( 2 . 71 ) 

i>{r, 6 ,t;e)(p/p)-P" = en? o ^ n ' ot ^T_P f (cos 0 )r-( f+1 ) + e( £ 2 ). ( 2 . 72 ) 

Order-of-magnitude analysis 

For this case the nondimensional total mechanical energy equation ( 2 . 43 ) reduces to 

■4 {K.E. + P.E.} = 0 . ( 2 . 73 ) 

at 
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The nondimensional kinetic and potential energies are, respectively, 


K.E. = / l|v| 2 dV + ~ [ hv | 2 dV 

(2.74) 

Jv m 2' pJv m T 

P.E. = / dS. 

(2.75) 


The averaging method requires the time-average of the kinetic and potential energies to be calculated to 
leading order. For this purpose, the potential flow approximations for the velocity field in each phase are 
used to calculate the leading order contributions to these quantities. Since the potential flow velocities in 
each phase are 0(e), the integrals over the material volumes in (2.74) may be expanded about a sphere of 
unit nondimensional radius to yield 

L 5 |v|i,i '' 

L I*** 

where v = eu and v = £u. 

Without small quantities in the integrand, the evaluation of the integral over the material surface in (2.75) 
is more complicated. Since the leading-order contribution from the kinetic energies is 0(e 2 ) the contributions 
from the potential energy are needed up to the same order. The assumed shape in (2.55) does not conserve 
volume to 0(e 2 ) and must be modified with an 0(e 2 ) correction in order to achieve this. In general, for a 
regular perturbation expansion of the nondimensional shape disturbance in the form 

r(0,t;e) = l + ef l (6,t)+e 2 f 2 (0,t) + O(e 3 ) (2.78) 


m -(ul 4- u 2 e )r 2 si n0drd6d(f 4- 0(s 3 ) 

m oo 1 

-(«£ +t)'r 2 sin OdrdOdip 4- 0(e 3 ) , 
2 


= £ 


= £ 


(2.76) 

(2.77) 


the volume of the shape can be shown to be 47t/ 3 -f 0(e 3 ) if f‘2(0,t) — -/?(#,£)• Using this modified shape 
the integral over the material surface may be expanded about, the unit sphere to yield 




1 + 2eh + 


•i 

£* 

T 



si nOdtpdO 4- 0(e 3 ) , 


(2.79) 


where fi(0,t) = }Pf(cos0). 

An alternate form for the time derivative of the nondimensional potential energy term, obtained with the 
use of the surface Reynolds Transport Theorem, is 


d 

dt 


j dS — J (V 5 • n)(n • v) dS . 


Here the unit normal and twice the mean nondimensional curvature are, respectively, 


n — e r — ee£-^-(0, t) 4- 0(e 2 ) 


V s n = 2 — e 


^(M)+cot0^(0,f)+2/,(M 


+ 0(e 2 ) • 


(2.80) 


(2.81) 

(2.82) 


By the divergence theorem and incompressibility, the constant part of the curvature term will not contribute 
to the surface integral on the right-hand side of (2.80). This allows the time derivative of the nondimensional 
potential energy to be expanded about a unit sphere in an alternate way 


d 

dt 




\d 2 h 

80' 2 


4- c°t + 2/i 


v r sin 0d0dip + 0(e 3 ) , 


(2.83) 


which does not require an 0(£ 2 ) correction to the shape. The use of the kinematic boundary condition (2.64) 
and an integration by parts shows that the time derivative of (2.79) is equivalent to (2.83). For consistency 
in the time-averaging, which is to be discussed next, (2.79) will be used throughout the chapter on theory 
to calculate the nondimensional surface potential energy. 
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Time-averaging 

Since the shape of the interface and the velocity solutions in each phase are complex functions, their real 
parts must be used in the integrands. The averaging method further requires only the time-average of the 
real part of the integrals making up the kinetic and potential energies. For example, the real part of the 
radial component of velocity in the drop is 

5ft{u r } = £P e (cosO)r e -' in to Ae iQtot + c.c.] + 0(e 2 ) . (2.84) 

Here c.c. refers to the complex conjugate. The nondirnensional time-average of the real part of the radial 
component of velocity squared is 

O, r2ir/Q(o 

<R{i; r } 2 ) = ^ / '«{<>}- r// (2.85) 

- 71 JO 

= e a P < 2 (cosg)r 2 ( < - 1) * n ” i4 * j +0(e 3 ). (2.86) 

Here | | 2 denotes the complex amplitude squared. This quantity may now be integrated over the sphere with 
unit radius to obtain part of the total contribution to the kinetic energy. 

\f I f (&{vr} 2 )r 2 sm ddrdedtp = |n m A| J + 0(£ 3 ) . (2.87) 

Similar operations, taking the real parts before squaring and time-averaging the integrands over one 
period, may be performed on all the terms in (2.74) and (2.75) to yield 


(K-E.) = 


(P.E.) 


e(e + + 1 ) 

, 2 V 

= 47r + e J 7r — 


’+!) + 


I^ro-dj 2 + 0(e 3 ) 


!)(^ + 2 ) | ^4 1 2 _)_ 0{e 3 ) . 


{2t + 1 ) 


(2.88) 

(2.89) 


Oscillator equation 

Following the averaging method, when a^(£) — }, the complex amplitudes, which represent time 

averages, may be replaced with terms involving a^(t) that would yield the same time average: 

|f^ce4| 2 — > 2 (df) 2 

|^| 2 -4 2 ( a ,) 2 . 

Rewriting (2.73) in terms of ae gives 


d_ 

dt 


e .,M±MM£l2( 4() > + 4, +e ^! 


Expanding and simplifying yields 

tU- l)(£ + l)(£ + 2) 

at + 


-!)(< + 2) 
{21 + 1 ) 


2 M 2 \ = 0 . 


at = O(e) . 


(2.90) 

(2.91) 


(2.92) 


[(£ + 1) + tp/p) 

Here the nondirnensional frequency for the shape oscillations is in agreement with Lamb’s result (2.59). 


(2.93) 


2.4.2 Case 1: “negligible” surfactant efFects 

This section discusses the viscous drop/medium system with no surfactant effects. The surface properties 
e*, «*, and p* are assumed to be “negligible”, or 0{a 3 ). All the surface dissipation and Marangoni terms 
on the right-hand side of the total mechanical energy equation (2.43) are neglected in this case. 

Due to the complicated forms for the equations and boundary conditions for the drop/medium system, 
this section calculates the uniformly valid velocity and pressure approximations to an accuracy of 0{sa). 
Unfortunately, this does not allow for the consistent expansion of the leading order nondirnensional total 
mechanical energy equation to 0(e 2 a 2 ). This will only affect the calculation of the resulting oscillation 
frequencies to 0(a 2 ), leaving the resulting calculation of the damping constant to 0(a 2 ) unaffected. 
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Inner equations 

Near the interface the radial coordinate may be replaced with r = R + y. Thus, y measures the distance away 
from the interface in the direction of the unit normal. Scaling the full Navier-Stokes equations in the drop 
using the inertial scales (2.42) with S = (p/ puj) 1 ^ 2 as the length scale for y introduces a small parameter 
a = S/R 7 the square root of the Reynolds number, a represents the ratio of the Stokes boundary layer 
thickness in the drop to the equilibrium radius. A similar Stokes boundary layer thickness 6 = (/t/pu>) 1/2 
and small parameter A = &/R is found in the medium phase. The nondimensional velocity (scaled with uiR) 
and pressure (scaled with pui 2 R 2 ) in the drop are written as regular perturbation expansions in a such that 

V(y,0,t,£,a) = £[V 0 (y,0,t) + aV 1 (y,e,t)] + O(£a 2 ) (2.94) 

P(y,0,t]£,a) = p eq +£[P 0 (y,e,t)+aP 1 (y,e,t)]+O{£a 2 ). (2.95) 

Similar regular expansions are assumed for the nondimensional velocity V and pressure P in the medium in 
terms of a. Here p eq (or p eq ) is the reference pressure in the drop (or medium) satisfying the nondimensional 
Young-Laplace equation p eq — peq(p/p) — 2. The size of the deformation is assumed to be small enough that. 
6 < a 2 and the fluid properties in each phase arc taken to be of the same order in magnitude. With the 
above expansions the inner equations in the drop phase take the following forms 



dVrO 

dy 

= 0 

(2.96) 


dP b 

dy 

= 0 

(2.97) 

dVg 0 dP 0 
dt dB 

O 

£ ^ 
i 

= 0, 

(2.98) 


with similar equations in the medium. These are the standard boundary layer equations for a flat surface 
driven by a pressure gradient supplied by the outer flow fields. The following boundary conditions hold at 
the interface 


Vro = Vro = Oho Ae^'Pt (cos 9) (2.99) 

Veo-Veo = 0 ( 2 . 100 ) 

dWo = 0 (2.101) 

dy a \/V dy 

p 0 - ?p 0 = (i- l){t + 2)Ae iilt ° l P t {cosB) . (2.102) 

P 


These represent, respectively, the kinematic, no-slip, tangential stress balance, and normal stress balance 
boundary conditions. The prescribed shape perturbation in (2.57) identically satisfies the normal stress 
balance boundary condition at this order. 

The equations for the 0(£a) inner variables take the following forms 


dVri 

dy 

m 

dy 

dVei m d 2 V<n 

dt dB dy 2 


-^r-cot0V e 0 -2Vro 
ov 

dVVo _ d 2 V r o 
dt dy 2 



+ 2 


dVeo 

dy 


(2.103) 

(2.104) 

(2.105) 


with similar equations in the medium. The 0(£a) variables are subject to the following boundary conditions 
at the interface 



Vn = Vrl 

= 0 

(2.106) 


Q. * 

Vo i Vo i 

a 

= o 

(2.107) 

dVm 

dy 

( fi\ l/i dV, n 

V p ) dy 

= 0 

(2.108) 
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Here the normal stress balance boundary condition is not enforced. Since the shape of the interface has been 
prescribed, not all of the boundary conditions can be satisfied to O(eo). As is customary in the analysis 
of drop deformations, see for example [30], the normal stress balance boundary condition is the degree of 
freedom that is lost. 

The boundary conditions far from the interface are replaced by matching conditions with outer velocity 
arid pressure solutions in each phase. For this problem, the leading order matching conditions state that 
the limit of the inner solution far from the interface must match the limit of the outer solution near the 
interface. For the higher order matching conditions a more precise definition is needed and the reader is 
referred to [54] for details. 

Outer equations 

Far from the interface the full Navier-Stokes equations are nondimensionalized using inertial scales with 
respect to the drop (2.42). If the nondimensional outer flow is irrotational and given by a scalar velocity 
potential (scaled with a;/? 2 ), where v = V0 in spherical axisymmetric coordinates, the viscous terms in the 
Navier-Stokes equations for both phases are identically zero. The pressure (scaled with pw 2 H 2 ) is given by 
the unsteady Bernoulli equation. Anticipating the matching conditions with the inner viscous flow fields, 
the outer potential flow fields are written as regular perturbation expansions in the parameter a 

0(r, 0,t]£,a) = e[<t>o{r,0,t) + a</>i(r, 0, t)] + 0{ea 2 ) (2.109) 

p(r,0,t;e,a) = p eq + e\po(r,0,t) 4- ap\ (r,Q, £)] + 0(ea 2 ) . (2.110) 

Similar expansions are assumed for and p in the medium in terms of a. 

The 0(e) and O(ea) outer equations are all governed by Laplace’s equation in spherical axisymmetric 
coordinates and the linearized unsteady Bernoulli equation for the pressure 

V 2 ^ = 0 

dfa 

Pi = -~m 

Similar equations for </> and p are used in the medium. The outer fields are subject to the boundary conditions 
that far from the interface the scalar velocity potential vanishes and the pressure tends to its reference value. 
The boundary conditions near the interface are replaced by matching conditions with inner velocity and 
pressure solutions in each phase. 

Time-periodic composite solutions 

An additive composite expansion [54] is used to obtain uniformly valid approximations for the time-periodic 
velocity and pressure fields in each phase. These take the following forms 

v tot = v(r, 9,t) + V (r, 0, t) — v m (r, 6, t) (2.112) 

p tot = p(r,0,t) + P(r,d,t) -p m (r,0,t) . (2.113) 

Here v tot and p tot represent the total composite fields and v m and p m represent the matched outer (or inner) 
fields in the overlap region. Similar expressions are used to construct the velocity and pressure in the medium. 
Solving the inner and outer equations in each phase, subject to their boundary and matching conditions, 
and constructing the uniformly valid solutions in outer variables yields the following approximations for the 
time-periodic velocity components and pressure. The approximate fields in the drop phase are given by 

i£ ot (r,0,t;e,a) = eififo.4e in,oi .Pf (costf) j/ -1 

- Q ^ 1 — (l + ^Ceolr*- 1 — EXP] 1 + 0(ea 2 ) , (2.114) 

(l + z)vWo J 

^ ot (r,0,£;e,a) = eiO <0 ^e in '°‘ t ^(cos6») |r f_1 + C eo (2 - ?-)EXP 
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!_ [(/ + lJCjor^ 1 - CfliEXP] 1 + 0(e a 2 ) , 

(1 + l) y/Ufo J 

p tot (r,0, t;e,a) — p eq = £n 2 to Ae iQtot ~P e (cosO)^r c 

+<1 (iTi)^ ( ' + 1)C “ r '} + 0( “ 2) ’ 

Similarly, the fields in the medium phase are given by 
fi‘ ot (r,M;e,d) 


v' 9 °'(r,6j;£,a) 


p tot (r,0j:£,a) - p eq 


Here EXP and EXP are functions that decay exponentially away from the interface in each respective phase. 
They are given explicitly by: 

EXP = exp [v/^ro ^ 1 ~ Z (2.120) 

L v2 ol 

EXP = exp -vAU (1 t i)(r : 1) l ■ (2121) 

v2 ol 

The real constants Coo , C$i , Coo , and C&\ are given in Appendix B. 

The above fields contain terms of two types: those that resemble potential flow or those that exponentially 
decay away from the interface. In the tangential component of the velocity, the exponentially decaying terms 
appear to leading order, or 0(e ), in the perturbed fields. In the normal component of the velocity, the 
exponentially decaying terms do not appear until second order, or O(sa). In the pressure, the exponential 
decay terms do not appear at all to this order in the expansion. 

The above velocity fields in each phase are used to approximate the terms in the energy equation. 

Order-of-magnitude analysis 

The direct substitution of the uniformly valid velocity approximations in the integrands of the nondimensional 
total mechanical energy equation leads to very complicated expressions. These expressions can be greatly 
simplified using an order-of-magnitude analysis to identify the most important terms. For this case the 
nondimensional total mechanical energy equation (2.43) reduces to 

4{K.E. + P.E.} = -a 2 {Bulk Diss.} . (2.122) 

dt 

The nondimensional energies and dissipation may be approximated using the uniformly valid velocity 
approximations from the matched asymptotic analysis. It is again convenient to use a new notation for the 


= £in io Ae iQeot P e (cos6) |r~ (f+2) 

1= =tCoo[r- (t+2) ~ EXP]\ + 0(ea 2 ) , (2.117) 

(1 + t) y/Qf 0 J 

= -einfflAe^'^y^fcosfl) jr (f+2 > + C eo ( 2 - r)EXP 

-d /t ^ - JL [^or- (<+2) +^ 1 ^]} + 0(ea 2 ), (2.118) 

(i + *)vn*o J 

= efl 2 o Ae' n ' ot ^-^yyP<(cos0) |r _(f+1) 

+O -/T * lC 9 0 r-< (+i q + O(g« 2 ). (2.119) 

(1 + i ) x/Sbo 


(2.115) 


(2.116) 
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(2.123) 


order-of-magnitudes for the nondimensional velocity in the drop phase: 

v tot = e[(u 0 + a Ul ) + (U 0 + aUt)] + 0(ea 2 ) , 


with a similar expression for the medium phase. Here the u’s refer to the potential flow terms in the 
uniformly valid velocity approximations with small 0(1) gradients in the radial direction. The volume 
integrals containing these terms will be significant over the nondimensional volume of 0(1). 

The U’s refer to the vortical flow terms in the uniformly valid velocity approximations with large 0{l/a) 
gradients in the radial direction that exponentially decay away from the interface. To distinguish these large 
gradients, the radial derivatives are recast in terms of the scaled variable y = (r - l)/a: 


dU 

dr 




1 d\J 

a dy 




(2.124) 


where dU/dy is an 0(1) quantity. The volume integrals containing the components of U are significant over 
the nondimensional volume of the boundary layer of O(a) near the interface in each phase. 

The following order-of-magnitude analysis is based on a substitution of these velocities into the integrands 
of the total mechanical energy equation. The integrals are written in spherical coordinates and expanded 
about a unit sphere. Only those contributions which are 0(e 2 a 2 ) and larger are kept. 

The nondimensional kinetic energy in the drop phase is 





+ u \ 0 ) + 2a(u r oU rl + ueouo i) | r 2 sin OdrdOdtp 
2(u r o| r==1 UrO + u 0o\ r =\ Ueo) + {UrO + Uio) 

* sin OdydOdif + 0(e 2 a 2 ) . 


(2.125) 

(2.126) 


A similar expansion for the nondimensional kinetic energy in the medium in terms of d. There the limits of 
integration of the potential flow velocity components are from r = 1 to r = oo, and the limits of integration 
for the vortical flow velocity components are from y = 0 to y = oo. Here the expansion has been truncated at 
0(s 2 a 2 ) because the uniformly valid velocity expansions accurate to 0(ea) do not contain the 0(ea 2 ) terms 
required to calculate the nondimensional kinetic energy to 0(£ 2 a 2 ). The level of algebra needed to construct 
the uniformly valid velocity to 0(ea 2 ) is enormous arid was not attempted. Neglecting the 0(£ 2 a 2 ) kinetic 
energy terms is shown below to affect only the calculation of the natural frequency of oscillations at 0(a 2 ). 
The calculation of the damping constant to 0(a 2 ) is unaffected. 

The nondimensional bulk viscous dissipation in the drop is 


L 


2(E tot : E tot ) dV = 


m i 

2 

. 


du 


rO 


dr 


+ “W 


1 f due o \ If , a \2 

( ~~qq~ “h ^ro j + ^ ( Wr0 — co ^ ^ U9 °) 


+ : 


dueo , 1 ( du r0 \] 2 1 2 

[d^ + r V“a0"" Uflo JJ Y 


sin 8drd0d(f 


n 2 7T r 7T /> 0 




/ 0 ^0 J -oo 


1 { dUgoY 1 (dugo du r0 

2 a 2 \ dy ) + a\ 9r + dO U °° 


dU, 


'00 


r—l 


dy 


a m 


sin Odydddip + 0(e 2 a ) . 


(2.127) 



Again, a similar expression is obtained for the bulk viscous dissipation in the medium in terms of d and its 
own limits of integration. Here the leading-order contribution to the nondimensional dissipation is 0(e 2 ot 1 ) 
and derives from the no-slip boundary condition at the interface. 


Time-averaging 

The averaging method proceeds by first multiplying the real parts of the velocity components making up the 
integrands, time-averaging these integrands over one period of oscillation, and finally integrating over the 
unit sphere. For the terms in the nondimensional total mechanical energy equation (2.122), this procedure 
yields 


(K.E.) = 


£(£ 4- 1)(2^ + 1) / 


+a 


(2£ + l) 2 


\fw 


u 

^fil 1 


+ 0(e 2 a 2 ) 


(P.E.) 
(Bulk Diss.) 


\[2 (s/VP + \fR>) 2 'fiho 

. .2 _ 1)(^ + 2 ) | 1 1 2 , / 

4 ’ + £ ’' (2f+l) Ml +< 

e 2 4 7T f 1 (2<+l) sfw 


£(£+1)\q 2\/2 (v/w+\/w) 

1 


o 


A 2 P 


2 + y/fip) 

+(i[(i + 2)p- {l- l)/5] 


2 (r - \)fip + 2£(e + 2) — 






A\ 2 + 0(e 2 a) . 


(2.128) 

(2.129) 


(2.130) 


The time-average of the total nondimensional potential energy has been calculated as in the base poten- 
tial flow using (2.79). Note again that the with the 0(ea) velocity solutions the time-average of the total 
nondimensional kinetic energy can only be calculated to 0(e 2 a). Since the time average of the total nondi- 
mensional viscous dissipation rate in the bulk is multiplied by a 2 in the nondimensional total mechanical 
energy equation, the kinetic energy term is the only quantity not known to at least 0(e 2 a 2 ). 


Oscillator equation 

Following the averaging method with «/>(£) = 3?{Ae l ^ f0 *}, the complex amplitudes in (2.128)-(2.130) are 
replaced with the quantities 

\tttoA\ 2 — > 2(a^) 2 (2.131) 

|A| 2 -> 2(^) 2 . (2.132) 


After simplifying, the total mechanical energy equation reduces to a damped harmonic oscillator equation 
of the form 

(1 -f- oAn) &i + (oBn 4- q 2 B^2) 4* fl 2 Q cii = 0 . (2.133) 

Here, 


An 
B n 


B^2 


(21 + l) 2 1 sfpp 1 

y/2 [(t + 1) + tp/p] (VPP + y/pp) V^O 

(21 4- l) 2 1 yfpp r^— 

\/2 [{£ 4 1)4 £p / p] {yfpp 4- > /pp) 

(21 + 1){2(1 2 - 1 )pp + 2£(£ + 2 )fi 2 p/p + p[(£ + 2)p - (£ - 1 )p]} 
[(£+!) + ip/ p](y/jlp+y/pp) 2 


(2.134) 

(2.135) 

(2.136) 


If the 0(a 2 ) added mass terms were calculated, the damped harmonic oscillator equation would take the 
form 

(1 4- aAn 4- q 2 A/ 2 ) 4- (orB^i 4- g 2 B^2) 4- cii = 0 . (2.137) 
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Assuming exponential behavior for the time-dependent oscillation amplitude of the form Uf (t) = e lUet , with 
Qf complex, the expressions for the frequency and damping constant can be found to be 


3{n*} 


n 


£0 


i An . 2 

( 3A 2 t 

B?, 


1 - a— — h a 

l 8 

1 

00 

JO , 
1 

o 

2 ). 


+ 0(c* 3 ) 


B 


— a- 


n 


+ ot‘ 


(AfiB^] — B«) 


+ 0(a 3 ). 


(2.138) 

(2.139) 


These show that contributes only to the expression for the frequency. Therefore, equation (2.133) can 
be expected to accurately predict frequencies to 0(a) and damping constants to 0(a 2 ). 

Interestingly, the results of (2.138) and (2.139) with A /2 — 0 agree with the 0[ol“) frequency and damping 
constant calculated by Marston [28], who solved the eigenvalue problem for the complex frequencies using a 
normal mode analysis accurate to 0{d 2 ). From this comparison it would appear that the 0(a 2 ) contribution 
to the added mass is zero. 


2.4.3 Case 2: “small” surfactant effects 

In this section, the viscous drop in vacuum is analyzed when surfactants are present. For this case, all surface 
terms on the right-hand side of the nondimensional total mechanical energy equation (2.43) are retained. 
The surface properties e*, k*, and p* are assumed to be “small 75 , or 0{a 2 ). 

This section presents the uniformly valid velocity and pressure approximations to 0(ea 2 ). This allows 
for the consistent expansion of the nondimensional total mechanical energy equation to 0(a 2 ), capturing 
both the frequency shift and damping times to this same order. 


Inner equations 

As in Case 1, the radial coordinate may be replaced with r — R + y in the inner region, and the full Navier- 
Stokes equations in the drop using the inertial scales (2.42) with the Stokes 5 boundary layer thickness as the 
length scale for y. The nondimensional velocity (scaled with ljR) and pressure (scaled with pcj 2 R 2 ) in the 
drop and the surfactant concentration at the interface (scaled with its equilibrium value T eq ) are written as 
regular perturbation expansions in a such that 

V(y,6,t\e,a) = £ [V 0 (y, 6, t) + aV, (y, 0, t) + a 2 V 2 (y,6,t)] + 0(ea 3 ) (2.140) 

P(y,Q,t\e,a) = p eq + e[P 0 (y, 0, t) + aPi(y,9, t) + a 2 P 2 (y,0,t)] + 0{ea 3 ) (2.141) 

r{0,t\e,a) = l+£[T 0 {e,t)+aT 1 (e,t)+a 2 T 2 {e,t)} + O{£a 3 ). (2.142) 


The size of the deformation is assumed to be small enough that £ « 2 - The inner equations in the drop 

phase take the following forms at O(e) 


dV B0 , dP 0 


dt 


+ 


00 


dT 0 w 1 fd 2 T n 

-df + Mo ~ pT.\'W + cot0 


U VrO 

dy 

dPo 

dy 

d 2 V e o 
dy 2 
dT o 
dO 


0 

0 

0 

0, 


(2.143) 

(2.144) 

(2.145) 

(2.146) 


where 


M o (0,t) 



+ cot 6 Veo 



y=0 


(2.147) 


These are the leading-order continuity, radial and tangential components of the linearized momentum equa- 
tions, and the surfactant transport equation for the perturbed fields. These fields are subject to the kinematic. 
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tangential stress balance, and normal stress balance boundary conditions at the interface 


Vro 

= iCl io Ae iQ, ° 1 P( (cos 9) 

(2.148) 

dVg o 

= 0 

(2.149) 

dy 



P b 

= (t- 1)(£ + 2)Ae ,n>ot P f (cos9) . 

(2.150) 


The prescribed shape in (2.55) identically satisfies the nondimensional normal stress balance boundary 
condition at this order. 

The equations for the 0(ea) inner variables take similar forms: 


dVoi dPi 

dt. oe 


dr, 


dVri 

dy 

m 

dy 

r e\_ 

dy 2 

or, 


d 2 V, 


d 2 v 


ot +M '-pb(^ L+co,e 8<> 


= -M 0 

dV r0 
dt 
dPo 

~ V ~d0 +2 dy 

= 0, 


rO 


dy 2 

dVg o 


where now 


M,(9,t) = 


dig, 

de 


+ cot 9 V 0 1 + 2V r , 


) 


y— o 


The 0(ea ) boundary conditions at the interface are 


Vn = 0 
dVg, dV r o 


dy 


_ , v e; dTo K ; dMo 
d9 60 a 2 d0 a 2 d9 


+ 




dN o 

de 


+ 2 cot 9 N 0 


P\ = 0 , 


where 


dV, 


90 

de 


N o (0,t) = 

Finally, the equations for the 0(ea 2 ) inner variables are 

dV r2 

dy 

dP2 

dy 

d 2 V 92 


COt 9 Vg 0 


) 


y=o 

-Mi + yM 0 

dV rl d 2 v ri 


dV 92 dP-2 
dt + d0 


dy 2 


dt dy 2 

dP\ , dVg, 2 dP 0 
V d9 + 2 dy V d0 


ar 2 i 

dt + 2 Pe s 


d 2 r. 
de 2 


where 


M 2 (9,t) 


(2.151) 

(2.152) 

(2.153) 

(2.154) 

(2.155) 

(2.156) 

(2.157) 

(2.158) 

(2.159) 

(2.160) 
(2.161) 


o dV 0o , dM 0 
2y dy + d9 

(2.162) 

Uco.^) = 0 

(2.163) 


(2.164) 


20 



The 0(ea 2 ) boundary conditions at the interface are again 


(2.165) 


V r2 = 0 

dV n _ dV Tl e; 5C < OM , 

Oy ~ 00 a 2 00 a 2 00 

+ 4f^- + 2cotflJV,V (2.166) 

a* \ OV ) 


where 


Ni (0,t) 


OVe i 
50 


cot OX 




9=0 


(2.167) 


The normal stress balance boundary condition is not enforced at 0{e a 2 ). This is similar to the loss of the 
normal stress balance boundary condition at 0(ea) in the viscous drop/medium system. The boundary 
conditions on the inner fields far from the interface are replaced by matching conditions with outer velocity 
and pressure solutions in each phase. 


Outer equations 

In the outer region far from the interface the full Navier-Stokes equations are nondimensionalized using 
inertial scales (2.42). Again, the viscous terms appear at G(a 2 ) but are identically zero for an outer potential 
flow. Anticipating the matching conditions with the inner viscous flow fields motivates writing the outer 
flow fields as regular perturbation expansions in the parameter a: 

4>(r,0,t;e,a) = s[4>o(r,0,t) + a<pi(r,8,t) + a 2 <p 2 {r,0,t)} + O(ea') (2.168) 

p(r,0,t]£,a) = Peg + e{po(r,0,t) + api(r,0,t) + a' 2 p 2 (r,0,t)\ + 0(ea 3 ) . (2.169) 

Here the arbitrary constant in the scalar velocity potential has been set to zero and p eq is the reference 
pressure in the drop satisfying the nondimensional Young-Laplace equation p eq = 2. 

With the above expansions, the G(e), 0(ea ), and 0{ea 2 ) outer equations are all governed by Laplace’s 
equation in spherical axisymmetric coordinates and the linearized unsteady Bernoulli equation for the pres- 
sure 

V 2 4>i = 0 

d(pi 

Pi = "IT 

These are subject to the boundary conditions that, far from the interface the scalar velocity potential vanishes 
and the pressure tends to its reference value. The boundary conditions near the interface are replaced by 
matching conditions with inner velocity and pressure solutions in the drop. 

Time-periodic composite solutions 

As in the Case 1, additive composite expansions are used to obtain uniformly valid approximations for 
the total time-periodic velocity field v tot and pressure field p tot in the drop. Solving the inner and outer 
equations in the drop, subject to their boundary and matching conditions, and constructing the uniformly 
valid solutions yields the following approximations for the time-periodic velocity components, pressure, and 
surfactant concentration written in outer variables: 

uj; 0t (r, 0, t; £, a) = si^l^Ae lUeQt Pt (cos 6) jr*” 1 

+ 0{ea z ) , (2.171) 

4 0t (r,M;e,a) - eiQmAe tn£ot j^{cosO) lr e ~ l + aCe\(2 - r)EXP 


, \/2 1 
(1 + i) VWo 


(i + l)Co\[r e ~ l - EXP] 


XP]j 


for i — 0, 1, 2 . 


(2.170) 
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(2.172) 


p tol {r,6,t;e,a) - p eq 


r (6,t;e,a) - 1 


-a 2 ,/* - JL [(l + 1 )C sl r ( - 1 - C 92 EXP] l + 0(ea 3 ) , 

(1 + t) y/Uto J 

£n%Ae int0, jP ( (cos6) |r f 

JL (< + l)C 91 r4 + 0(ea 3 ) (2.173) 

(1 + 0 v/Hfo J 

eGe m ‘ ot |l + aCgi + a 2 C 92 j Pe(cos9) + 0(ea 3 ) . (2.174) 


Here EXP is an exponentially decaying function given by 


EXP = exp 



(1 + t) (r - 1) 
y/2 a 


(2.175) 


The complex constants Co i, C021 G, C 5 i, and C g 2 are given in the Appendix B. 

Again, the above fields contain terms of two types: those that resemble potential flow or those that 
exponentially decay away from the interface. In the tangential component of the velocity, the exponential 
decay terms appear at 0(eol) in the perturbed fields, which is one order higher in a that in the viscous 
drop/medium system of Case 1. In the normal component of the velocity, the exponentially decaying terms 
do not appear until 0(ea 2 ). In the pressure, the exponentially decaying terms do not appear at all to this 
order. 


Order-of-magnitude analysis 

As before, the direct substitution of the uniformly valid velocity and surfactant concentration approxima- 
tions in the integrands of the nondimensional total mechanical energy equation leads to very complicated 
expressions that can be greatly simplified using an order-of-magnitude analysis. Recall that for this case the 
nondimensional total mechanical energy equation is 


s' KE+PE) 


The Marangoni term can be written as 

{Marangoni} 


— « 2 {Bulk Diss.} 
— p*{Surf. Diss.} 
+e* {Marangoni} . 


-4* {Mar. E.} 
dt 1 ; 

__L_{Mar. Diss.} 

rGj 

— {Remainder} . 


(2.176) 


(2.177) 


Equations (2.45)-(2.53), with the fluid properties in the medium set to zero, define the quantities in curly 
brackets. 

The nondimensional energies, dissipation rates, and Marangoni terms may be approximated using the 
uniformly valid velocity approximations from the matched asymptotic analysis. As before, it is convenient 
to use a new notation for the order-of-magnitude of the velocity in the drop: 

v' 0 ' = e [(uo + aui + a 2 u 2 ) + (U 0 + aU, + a 2 U 2 )] + Q(ea 3 ) . (2.178) 

For this case, the only nonzero components in the expansion (2.178) are u r0 , u r2 , ug 0 , ug 2 , U r2 , Vei, and 
U e 2 . Here the u’s refer to the terms resembling potential flow with small 0(1) gradients. The volume 
integrals containing these terms will be significant over the nondimensional volume of 0(1). The U’s refer 
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to the vortical flow terms with large 0{\ja) gradients in the radial direction that decay exponentially away 
from the interface. The volume integrals containing these terms will be significant over the nondimensional 
volume of the boundary layer of 0{a) near the interface. 

Following the order-of-magnitude analysis of Case 1, keeping only those contributions which are 0{e 2 a 2 ) 
and larger, the nondimensional kinetic energy in the drop is 



dV = 



4- u 2 e 0 ) + 2a 2 (u r0 u r2 + U 0 QU 02 ) | r 2 sin 0drd8d<p 
2a u eo \ r=l Uo\ \ sin OdydOdtp + 0(s 2 a 3 ) . 


(2.179) 


Unlike the previous drop/medium case, all the nondimensional kinetic energy terms to 0(e 2 a“) may be 
calculated because the uniformly valid velocity profiles have been determined to O(ea'). This will allow for 
the calculation of the frequency shift to 0(a 2 ). 

The bulk viscous dissipation rate in the drop is 


L 2(E 


tot . E totx d y _ 


2 IT PIT /•] 


( du r 


0 Jo Jo 


\ dr ) + r 2 


1 ( duQQ 


dO 


+ W r0 + — (UrO ~ COtOu 0 o) 


1 

+ 2 


duoo , 1 ( du r o 

H 1 ”7^ u 90 


)]•} 


r 2 sin OdrdOd^f 


dr r \ dO 
+0(e 2 a ) . 

The surface viscous dissipation rate arising from the surface shear viscosity is 


p p2tt pit 

/ 2(E* s ot : E* ot ) dS ~ e 2 / / n 2 sin BdBdup + 0(e 2 a 

J Jo Jo 


where 


n 0 (9,t) 


_ ( dug 0 

“ V 99 


cot 0 uo 0 


The surface viscous dissipation rate arising from the surface dilatational viscosity is similarly: 


p p2tt ptt 

/ (V s • v tot ) 2 dS = e 2 / ml sin OdOdip -F 0(£ 2 a ) , 

JSm J° J° 


where 


f r\ 

m 0 (6, t) = ( + cot 9 ugo + 2u r0 


89 


r= 1 


The nondimensional Marangoni term is 

[ (r - 1 )(V s ■ v tot ) dS - e 2 f f r oirio sin QdQdip + 0(e 2 a) . 
Js m Jo Jo 

The Marangoni stored energy term is 

f hr - l) 2 dS = e 2 f f ho sin OdOdip + 0(e 2 a) . 

Js,„ 2 Jo Jo 2 


(2.180) 


(2.181) 


(2.182) 


(2.183) 


(2.184) 


(2.185) 


(2.186) 
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(2.187) 


Finally, the leading-order Marangoni dissipation rate is given by 

j |V.(r-l)| 2 dS = e 2 J^ ^ (^) sin 0d6dip + 0(e 2 a) . 

The expressions in (2.186) and (2.187) need to be calculated only to 0(£ 2 ) since they couple with the 
nondimensional total mechanical energy equation through the Marangoni term, which is multiplied by e*, 
an 0(a 2 ) parameter for this case. 

The integrand of the Remainder term is 0(e 3 ) and hence neglected. 


Time-averaging 


Using the time-averaging procedure as outlined in Case 1, the time-average of the terms in the total me- 
chanical energy equation (2.176) and the Marangoni expression (2.177) are given by 


(K.E.) 

(P.E.) 

(Bulk Diss.) 
(Surf. Diss.) 

(Marangoni) 
(M.E.) 
(Mar. Diss.) 


£ 2 7T 


1 


£(2^+1) 


|fWl| 2 + 0{eW) 


4 ir + e 2 J e 1 )( £ + 2 )i ,112 , 


“47T 


(2i + 1) 
('-D.O ,02 


-lAf + CHe 3 ) 


l 


ftf 0 .4| 2 + 0(e J a) 


= e 27 r 


2o_ U 


t(2t + 1) 
+0(£ 2 a) 


(£ + l)(f + 2) + t(t — 1)— j 

f l S 


|« f0 .4| 2 


= €7T 


(t- i) 


{in? 0 AG* -(- c.c.) + (D(£ 2 ct) 


(21+1) 

^WTn' G ' 7 + 0 ^ 

£22 'i^T' G ' s + 0 ' fi “>- 


(2.188) 

(2.189) 

(2.190) 

(2.191) 

(2.192) 

(2.193) 

(2.194) 


Coupled oscillator equations 

Following the averaging method with a e (t) = 'R{Ae iniot } and g e (t) = 9?{GV fi '° r }, the complex amplitudes 
in (2.188)-(2.194) are replaced with the quantities 


|ftro.4| 2 — t 2(d()~ (2.195) 

| 2 l| 2 -> 2 (a e ) 2 (2.196) 

+ c.c.) — > 2(g(d{) (2.197) 

|G| 2 -4 2(<?<) 2 . (2.198) 


These would give the same time average if a/(t) and g e (t) were time-periodic with a nondimensional period 
2n/Cl{o. Simplifying the nondimensional total mechanical energy equation yields the damped harmonic 
oscillator equation for the nondimensional amplitude of the shape oscillations, 

a t + d e [a 2 2(£-l)(2e + l) + K’ a e(e-l) 2 + g* s (i-l)(e + l)(£ + 2)} 

+ at[l{l -m + 2)] = - fl< [e:f(f-l)]. (2-199) 

This oscillator equation is coupled to the simplified Marangoni expression: 

9t + y^9t = {t- Do*. (2-200) 

1 c s 

Setting a 2 ,e*,tt*, and /i* to zero reduces the coupled set of equations to a simple harmonic oscillator 
equation for a^(£), the time-dependent amplitude for the shape oscillations of a drop in vacuum. The 
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oscillation frequency is in agreement with the results of Lamb [23] when the properties of the medium are 
negligible compared to those of the drop. 

By setting e* — 0 the equations decouple and the oscillator equation for at(t) contains damping contri- 
butions from a 2 ,K*, and //*. 

For large surface Peclet number Pe s , the time-dependence for the surfactant concentration may be seen 
to be proportional to the time-dependence for the amplitude of the shape oscillations. Therefore the term 
on the right-hand side of (2.199), which arises from the Marangoni effect, acts to “stiffen 5 ' the system and 
increase the natural frequency of oscillation. 

For small Pe s , equation (2.200) shows that gt(t) is proportional to a f \t). The term on the right-hand side 
of (2.199) thus gives rise to an additional damping mechanism from the irreversible diffusion of surfactant 
molecules. 

For a surface Peclet number of 0(1), the coupled system (2.199) and (2.200) needs to be solved simulta- 
neously to describe the oscillations. 


2.4.4 Case 3: “medium” surfactant effects 

In this section the presence of the surrounding medium is neglected and the viscous drop in vacuum system 
when surfactants are present is analyzed. For this case, all terms on the right-hand side of the nondimensional 
total mechanical energy equation (2.43) are retained. The surface properties e*, k*, and /i* are assumed to 
be “medium”, or 0(a). 

Due to the larger influence of the surfactant monolayer, the uniformly valid velocity and pressure ap- 
proximations are calculated to an accuracy of 0{ea). This allows for a consistent expansion of the total 
mechanical energy equation to 0(a), capturing both the frequency shift and damping times to this same 
order. 


Inner equations 

As before, the inner fields are written as regular perturbation expansions in a in the scaled variable y: 

V(y,6,t;e,a) = e{y 0 (y,$,t) +aV 1 (y,e,t)] + 0(ea?) (2.201) 

P{y,0,t-,£,a) = Peg +e[Po(y,&,t) + aPi{y,9,t)] + &(ea 2 ) (2.202) 

r(M;£,<*) = 1 +£[r o (0,t) d-aT^.t)] + 0{ea 2 ) (2.203) 

With the above expansions, the 0(e ) and 0(eoc) inner equations take the same forms as in Case 2 equations 
(2.143)-(2.146) and (2.151)-(2.154), respectively. The kinematic, tangential stress balance, and normal stress 
balance boundary conditions at 0(e) are 


V r o 

dVeo 

dy 

Po 


m eo Ae incot P e {cos6) 


e* s dT 0 k* s dMp nl (dNp 

a 89 a 89 a \ 89 

{e-l){e + 2)Ae iQcot P e (cos9). 


+ 2cot ,9 


Ao) 


(2.204) 

(2.205) 

(2.206) 


N 0 given by (2.159). The prescribed shape (2.55) identically satisfies the normal stress balance boundary 
condition. 

The corresponding boundary conditions at 0{ea) are 

V r i = 0 (2.207) 

8V 61 8Vr o , e ; sr, t K ; 8M , 

dy ~ 89 + a 89 + a 89 

+ ~ (^ + 2cot9N 0 ). (2-208) 

a \ o0 ) 

The normal stress balance boundary condition cannot be satisfied at this order. 
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Time-periodic composite solutions 

The outer equations are identical to those from Case 1. Thus, the composite expansions may be obtained 
as before and the composite expansions are given by 


vl ot (r,0,t;e,a) = 


eifl eo Ae ineot Pe(cos6) |r <_1 

V2 1 


+a 


vg ot (r,0,t\£,a) = 


(1 + i) \ZTTe o 

n 

d e 


{t + VCeoV 1 - 1 + EXP] } + 0(ea 2 ) , 


(2.209) 


ein eo Ae tn ' ai ]^§-( cos 6) jr*" 1 + C eo (2 - r)EXP 


-a ^ ^- [(t + 1)CW/-' - EXP] \ + 0(ea 2 ) , 


(2.210) 


en 2 eo Ae in ‘ o1 -P f (cos 6) \r e 


p iot (r,9,t;e,a) — p eq = 


(1 + i) s/iieo 

s 

r{0,t\e,a)-l = £Ge in ‘ ot ( 1 + aC sl 1 A(cos0) + C(ea 2 ) , 


-aTT^T 7 -i(f + 1 )C eo r e 1 + 0(ea 2 ) 


( 2 . 211 ) 

( 2 . 212 ) 


EXP is the function in (2.120) that decays exponentially away from the interface in the drop. The complex 
constants Coo, Co 1 , G, and C 9 \ are given in the Appendix B. 

The fields in (2.209)-(2.211) take the same form as the the fields in the drop phase (2.114) (2.116) for 
the drop/medium system considered in Case 1. In Case 1, however, the presence of the 0{e) vortical field 
in the tangential component is traced to the viscous forces between the fluids in each phase in satisfying the 
no-slip boundary condition at the interface. In this case, on the other hand, its presence comes from the 
forces exerted between the fluid in the drop and the surfactant monolayer. 


Order-of-magnitude analysis 

The nondimensional total mechanical energy equation and Marangoni expression are identical to Case 2 
equations (2.176) and (2.177). 

The energies, dissipation rates and Marangoni terms may be approximated using the uniformly valid 
velocity approximations from the matched asymptotic analysis. It is again convenient to use a slightly 
different notation for the order-of-magnitude of the velocity in the drop: 

v tot = e[(u 0 + qui) + (Uo + qUi)] + 0(ea 2 ) . (2.213) 

For this case, only the U r o component is zero. As before, the u’s refer to the terms resembling potential 
flow with small (9(1) gradients. The volume integrals containing these terms will be significant over the 
nondimensional volume of 0(1). The U’s refer to the vortical flow terms with large 0(1 /a) gradients in 
the radial direction that decay exponentially away from the interface. The volume integrals containing these 
terms will be significant over the nondimensional volume of the boundary layer of O(a) near the interface. 

Substituting the composite expansion into the expression for the kinetic energy and expanding to 0(e 2 a) 
yields: 



dV = 


e 


2 



(u 2 ro + u 2 0o ) + 2a(u r0 Url + UQoUo i) 


r 2 sin OdrdOdip 
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nr i 

JO Jo J — OG L 


+e 2 a 


2 '2 x 


2 w 0o | r=1 Ueo + C/| 0 f sin OdydOdtp + 0(6* a 


(2.214) 


Similarly, the bulk viscous dissipation rate in the drop is given by: 

[ 2(E tot : E tot ) dV = 

J v m 


m o 
- 


£ 2 Q 


1 /5L r ( 


2q 2 \ dy 


Oo 


sin QdydOdtp + 0(e 2 ) . 


The surface viscous dissipation rate arising from the surface shear viscosity: 

p p2lX p 7T 

/ 2(E^ ot : E^ ot ) dS = e 2 / (^o + No) 2 sin0d0dy? + 0(£ 2 a) . 

*/o io 

Here, n 0 and Nq are given by: 


n o (0,t) = 


du 


eo 


do 


— cot Ouqq 


N O (0,t) = (^-COtOUgo 


r=l 


r=l 


The surface viscous dissipation rate arising from the surface dilatational viscosity is similarly: 

p p2i\ pi r 

/ ( V s • v tot ) 2 dS = £ 2 j / (m 0 + Mo) 2 sill OdOdif + 0(e 2 a) , 

Js m J o Jo 


where 


m 0 (6,t) = + cot^o + 2u r0 


M o (0,t) 


f^T +cot OUm + Wrt 
\ 90 j 


r~ 1 


r = l 


The nondimensional Marangoni term is given by 

p2l[ p 7T 


P pzi\ pi\ 

/ (r - 1)( V, • v tot ) dS = e 2 / r 0 (TO 0 + Mo) sin OdBdip + 0(e 2 a) . 

Js„ Jo Jo 


The stored Marangoni energy is 


f -(r - l) 2 dS = e 2 [ [ ^Tq sinOdOdip + 0(e 2 a ) . 
J s m 2 Jo Jo * 


Finally, the Marangoni dissipation rate is 


p P 2tt p tt / ap \ ^ 

/ |V,(r - 1)| 2 dS = £ 2 / ( ) sin OdOdtp + (9(e" 

Js m Jo Jo \ do / 


a). 


(2.215) 

(2.216) 

(2.217) 

(2.218) 

(2.219) 

( 2 . 220 ) 
( 2 . 221 ) 

( 2 . 222 ) 

(2.223) 

(2.224) 


The integrand of the Remainder term is 0(e 3 ) and again neglected. 


Time-averaging 

As before, the order-of-magnitude estimates for the integrals in the energy equation are time-averaged over 
one period of oscillation and integrated over the unit sphere. For the terms in the nondimensional total 
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mechanical energy equation (2.176) and the Marangoni expression (2.177), this procedure yields 


(K.E.) = e 2 7T 

+ 


1 


*(2*+ 1) 

2t 


1 + 


a 


+ 1 ) 


y/n^l{2l+ 1) 


ICVol 2 


Vi 


(Cooe 1,r/ ^ 4 + c.c.) 


(P.E.) = 47T + £ J 7r 


(2*+l) 

2 je- m + 2) 




|flfo^4| + 0(e a ) 


( 21 + 1 ) 


\A \ 2 + 0(e 3 ) 


(Bulk Diss.) = -47r^-^|C 90 | 2 |^o^r 2 +0(e 2 ) 

a i 


{Surf. Diss.) = £ 2 27r 

+ 


(P-i) ne + 2) 
(2t + i) \ e 
(<-i) 


[l + (Coo + Coo) 4- ICpol" 


(*+D 


+ (Coo + Coo) + l^l^ol 2 


-T 

! l s 1 


(2.225) 

(2.226) 
(2.227) 


+0(s 2 a) 

(2.228) 

(Marangoni) = e 2 n + ^ [(£ l)(ifi* 0 .AG + c.c.) 


+ {£ + l)(in^o^o AG T c.c.)] + 0(s 2 a) 

(2.229) 

(Mar. E.) = e^—^Gf + O^a) 

(2.230) 

(Mar. Diss.) = £ 2 2 t r + j |G| 2 + 0(e 2 a) . 

(2.231) 

Coupled oscillator equations 


Following the averaging method with at(t) = and gi(t) = 

in (2.225)— (2.231) are replaced with the quantities 

fo( }, the complex amplitudes 

1 0 ^4 1 2 — > 2(a^) 2 

(2.232) 

cT 

t 

C<i 

(2.233) 

(iQto-AG + c.c.) — t 2 Qtat 

(2.234) 

|G| 2 -> 2(^) 2 

(2.235) 


which would give the same time average if c if(t ) and gt(t) were time-periodic with a nondimensional period 

27r/n^ 0 - 

Due to the presence of the complex constant Coo , it is not clear how to replace the time-averaged quantity 
(iVlioCooAG + c.c.) in the Marangoni term with terms involving d/ and gi. A result may still be extracted 
from the above analysis if the Gibbs elasticity is set to zero. The limit when e s = 0 corresponds to the case 
when there is so much insoluble surfactant on the interface that local changes in shape do not change the 
surfactant concentration appreciably. For this situation the surfactant concentration and surface tension are 
approximately constant over the interface. The Marangoni term that couples the oscillator equation for the 
time-dependent amplitude of the shape deformation to the surfactant transport equation then vanishes. The 
resulting decoupled oscillator equation for at(t) is given by 


{ 


df \ 1 + 


+ at < a 


is/flto 


a (*+l) 
1(21+1) 
(*+l) 


\Ceo\ 2 


+ 


21 


s/2 (2£ + 1 ) 


(Cg oe 


-i tt /4 


+ C.C. 


} 


s/2 


-|C«o I 2 


+/x*(£ — 1)(^ + 1)(^ + 2) [l + (Cqo + c.c.) + |C r <9o| 2 ] 

+K*J [(£ - l ) 2 + (t- l)(e + l)(Cm + c.c.) + (t + l) 2 |C ro | 2 ] J 
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(2.236) 


+ atW- l)(* + 2)] 

= 0(a 2 ) . 

Here the complex constant Cgo is independent of e* since it has been set to zero. Although the expressions 
for the 0(a) added mass and damping terms may be found by evaluating the real terms, they are extremely 
complicated and are not written explicitly here. 

If the surface viscosities are zero, the uncoupled oscillator equation reduces to 

a t + [e(e-l)(f. + 2)}a e = 0(a 2 ) (2.237) 

This limit shows that without surface viscosity, the 0{a) added mass and damping terms are shifted to 0(a 2 ), 
which is neglected in this leading-order analysis. Without the surface viscosities or Marangoni effects, the 
viscous bulk fluid in the drop has nothing to “push against” to produce O(o) dissipation in the boundary 
layers. 

If the bulk viscosity of the drop is neglected, the oscillator equation simplifies to 

. f Al(t-W + 2 1 

a ' + 0< \[(/-i)(/ + 2K + ^ + !)«;]/ 

+ a* [1(1 - 1)(* + 2)] = 0(a 2 ) . (2.238) 

Here a damping term remains from the dissipation within the interface due to the surface viscosities. Note 
that in the limit when one of the surface viscosities is larger than the other, it is the smaller of the two 
surface viscosities which contributes the most to the damping. The natural frequency of oscillation remains 
unchanged from the Lamb frequency to 0(a), but the presence of the 0(a) damping term in this limit causes 
a frequency shift at 0(a 2 ), which has been neglected. 


2.4.5 Case 4: “large” surfactant effects 

In this section, the leading-order effects of viscosity and surfactants on the drop/medium system are analyzed 
when the surface properties ej, fi* s , and k* are assumed to be “large”, or 0(1). 

For this case the viscous dissipation terms in the bulk phases are negligible compared to the surface 
viscous dissipation and Marangoni effects at the interface. For this reason the uniformly valid velocity 
profiles in each phase and the analysis of the energy equation using the averaging method is not needed. The 
dynamics of the drop is determined by the potential flow, but potential flow cannot satisfy the tangential 
stress balance boundary condition. To treat this, V$ y the tangential velocity of the interface is introduced. 
Vq is distinct from the tangential velocity of the fluid in each phase evaluated at the interface. 

The assumed solutions and leading-order field equations are thus of the form: 


0(r, 0 , £; e, a) 

- e<t> 0 (r,d, t) + 0(ea) 

(2.239) 

p(r,Q,t',e,a) 

= p e g + ep 0 (r,6,t) + 0(ea) 

(2.240) 

4>(r,0,t;e,a) 

= e<f>o(r, 0, t.) + O(ea) 

(2.241) 

p(r , 6, £;e, a) 

= Peq + epo(r, 6, t ) + 0{ea) 

(2.242) 

T(0,£;£,a) 

= 1+er O (0,t) + 0(ea) 

(2.243) 

Vo(6,t\e,a) 

= eVga{9, t) + O{eo) 

(2.244) 


5To 

dt 


V 2 <£o = 0 , po = 

for r < 1 

(2.245) 

o - p <9</>o 

V 2 0 O = 0 , * = -*-£ 

for r > 1 

(2.246) 

1 ( 92T °^ 

= Pe, U 2 +C ° t# W j 

at r = 1 . 

(2.247) 
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These are subject to the O(s) kinematic, tangential stress balance, and normal stress balance boundary 
conditions at r = 1 , respectively, 


dtpo _ chpo _ dj_ 

dr dr dt 

= e; ar 0 k* s dMp ^ 

a dO a d6 

»-** = -{U 

In the above expressions Mq and Nq are given by 

M 0 (0, t) = + coteVeo + 2u r0 ) 

N 0 (6,t) = (j^--cot6V 0o y 


4 (^ + 2«#*) 


r— 1 


(2.248) 

(2.249) 

(2.250) 

(2.251) 

(2.252) 


These contain the tangential velocity of the interface V^o and the normal velocity velocity of the fluid 
evaluated at the interface i> r o = d(j)/dr. 

The surface deformation f(0,t ), surfactant concentration To and the tangential velocity of the 
interface V$o are assumed to have the forms 


f(0,t) 

= at(t)Pi(cos8) 

(2.253) 

r 0 (o,t) 

= g e {t)P e {cos6) 

(2.254) 

Voo 

dPi 

= b e (t)--J-(cos6). 

(2.255) 


Here a?, gt, and be are the time-dependent amplitude perturbations for a pure mode. The kinematic boundary 
condition (2.248) allows a similar modal expansion to be obtained for the velocity potentials 0 o and 0 o with 
amplitudes proportional to a?. 


Coupled oscillator equations 

Substituting the modal expansions (2.253)-(2.255) into the normal stress balance boundary condition (2.250), 
the tangential stress balance boundary condition (2.249), and the surfactant transport equation (2.247) leads 
to the three coupled equations: 



P 1 
P (*+ 1) 


at 


0 


gt — £(£ T 1)6/ T 2 o>i 


-0 1 - l)(t + 2 )ae + 2 e* s gt - <[4 d e + 2£(£ + 1 )6,] 

-e; 9 t + «;[2d/ - ea + 1 )b t ] + ti(e- m + 2)6/ 
t(i + 1) 

r>_ ■ 


Using equation (2.257) 


to eliminate bi(t) leaves the two coupled equations: 


dt[{£ + 1) + tp/ p] 


■ f An'Mt- W+W + 2) 1 
a T[«^ + i) + M:(^-i)(^ + 2)]j 

+ at [t(t — 1)(£ + 1)(^ + 2)] 

f 2e;fi;t(t-i)(t + i)(e + 2) ] 
5f \ [«;/(/+!)+/*:(< -!)(< + 2)]/ 
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f 1(1+1) e;i(/+l) 
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(2.257) 

(2.258) 


(2.259) 


(2.260) 
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Equations (2.259) and (2.260) describe the leading-order dynamics of the system. They neglect all effects at 

0(a)- 

If the surface shear viscosity p* s is set to zero, the oscillator equation for the amplitude of the shape 
oscillations (2.259) decouples from the surfactant transport equation and reduces to 


di + a? 


1 ( 1 - l)(l + l)(l + 2) 

[(e + i) + ep/p\ 


= 0(a) ■ 


(2.261) 


Thus, even when the surface dilatationai viscosity is nonzero, or the surface tension depends on surfactant 
concentration, the leading-order dynamics for the shape oscillations of the system are, in the absence of 
surface shear viscosity, the same as those for an inviscid drop/medium system. 

If the Gibbs elasticity is zero, equations (2.259) and (2.260) decouple. The oscillator equation (2.259) 
thus reduces to 


4* 1) + tp/ p\ 


. r akwm-w + w + 2 ) \ 

af \[K^ + i) + p;(<?-i)(* + 2)]/ 

+ a e (£(£-m + m + 2)] = 0(a) 


(2.262) 


Equation (2.262) is valid for a drop/medium system and generalizes the result (2.238) obtained in Case 3, 
which neglected the bulk viscosities and Gibbs elasticity for a drop in vacuum. 
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Chapter 3 


Numerics 


3.1 Governing equations and boundary conditions 

Here the potential flow solutions for moderate-amplitude shape oscillations of a drop in vacuum without 
gravity are considered. The inviscid drop has an unperturbed equilibrium radius R , a density p, and a clean 
interface with constant surface tension a Q . 

As before, the Reynolds number associated with the shape oscillations is assumed large enough that 
deviations from inviscid flow are present only in thin boundary layers near the surface. This implies that 
the fluid motion in the bulk of the drop is well approximated as inviscid. For the moderate deformations 
considered here, the flow r is further assumed to be irrotational. The governing equation is then Laplace’s 
equation for a scalar velocity potential (p. 

Because of their ability to track the position of the interface, boundary integral methods are particularly 
well suited to this free-surface problem. There are two common boundary integral formulations for Laplace’s 
equation: the single-layer potential representation and the double-layer potential representation. The latter 
has been shown [42] to be a more accurate and efficient formulation for the numerical simulation of the 
shape oscillations of two-dimensional drops. This section discusses the extension of the double-layer potential 
boundary integral formulation to the case of axisymmetric shape oscillations. 

3.1.1 Regularized boundary integral representation 

The free-space Green’s function for Laplace’s equation in three dimensions is 

G(xi,x) = -T“ T — 7 with V 2 G(x;, x) = <$(x - Xj) . (3.1) 

4tt |x — Xj| 

x is the field point position vector and x t is the source point position vector. 

For a source point x* located on the boundary of the domain, the double-layer potential boundary integral 
representation [35] 

[ n(x) • VG(xi,x) g(x) dS(x) + ^(x*) = ^(xj) , (3.2) 

JS(x) 1 

expresses the harmonic scalar velocity potential (p in terms of a surface distribution of dipole strengths q 
oriented parallel to the outward pointing unit normal vector n. Specifying the scalar velocity potential (p 
reduces the double-layer representation (3.2) to a Fredholm integral equation of the second kind for the 
dipole strengths q. The dipole density distribution may be used to find a vector velocity potential through 
the integral relation [35] 

•0(xi) = - [ n(x) x VG(x;, x) q(x) dS(x ) . (3.3) 

•'S(x) 

The velocity of the interface can be expressed by the two distinct representations: v = Vcp = Vx^. 
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Figure 3.1: a) The local coordinate system in the meridian plane with variable parameterized with arclength 
s. b) The local coordinate system in a plane perpendicular to the axis of symmetry. 


The identities [42] 


f n • VG(x n x) dS(x) 

JS(x) 

f n x V(7(xj, x) dS(x) 

Js(x ) 


1 

2 

0 , 


(3.4) 

(3.5) 


may be used to express the integral representations and relations in their regularized forms 



* VG(x,,x) [<?(x) - q{x t )]dS(x) 


+ <?(*;) 


0(x z ) 


(3.6) 


ip(xi) = - [ h(x) x VG(xi,x)[q{x) - q(xi)]dS(x) . (3.7) 

JS{x) 

The integrands in (3.6) and (3.7) are non-singular everywhere in the integration domain S(x). The regular- 
ization of the integral representations and relations, along with the choice of a convenient local coordinate 
system, to be discussed next, greatly simplifies the numerical problem. 

The boundary integral formulation, which reduces the three-dimensional problem to a two-dimensional 
one, may be further simplified with the assumption of an axisymmetric geometry, which reduces the two- 
dimensional integrals over the surface to one-dimensional line-integrals over the contour of the domain in a 
meridian plane. For convenience, a local coordinate system is introduced as illustrated in Figure (3.1). The 
position vector x(s,<^) = e r (</?)r(s) +e z z(s) and all other variables are parameterized in terms of arclength 
s, where ds 2 = d x* dx for a fixed meridian angle <p. Note that r represents the distance from the z-axis and 
not that from the origin. Explicit relations for this local coordinate system are given in Appendix C. For the 
axisymmetric case, all variables except the unit vectors e r and e^, are independent of c p. The infinitesimal 
surface area is dS — r dsd^p and the (^-dependence may be integrated analytically to give the parameterized 
forms ^ 

f h(s) ■ VG(.Si, s) r(s) [ q{s ) - <?(sj)] ds + q(s ,) = <£(«,•) (3.8) 

Jo 

ip{ Si ) = - f * [n(s) x VG(si, s) r(s)] [q{s) - q(s l )} ds , (3.9) 

Jo 

L is the total arclength between the “poles” of the axisymmetric drop, is evaluated in the meridian plane 
which contains x*. The vector velocity potential is ip = e^ip. The axisymmetric kernels in (3.8) and (3.9) 
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can be expressed in terms of the complete elliptic integrals of the first and second kinds. Explicit forms for 
these weakly singular kernels and their asymptotic behavior as s -> Si may be found in the Appendix D. 
Unlike the regularized integrand in (3.8), which is zero at s = s t , analysis of the integrand in (3.9) reveals 
that as s -> Si 

e v ■ [h(s) x VG(sj,s)r(s)] [q(s) - g(si)] ~ (Sz) + ~ Si ) * (3- 10 ) 

Since v = V</> = V x ip the velocity of the interface v = hv n + sv s has components that may be 


calculated from derivatives in arclength only 



X 

> 

II 

£ 

= 2<*> + $*<•> 

(3.11) 

v s (s) - s • V</> 

- £<•>■ 

(3.12) 

3.1.2 Boundary conditions 




The potential flow solution to the integral representation in (3.8) must satisfy the kinematic and normal 
stress balance boundary conditions 


= n v n (3.13) 

= <T efl (C,+C v ). (3.14) 

Here, n(s) = -e r z'(s) + e 2 r'(s) is the outward pointing unit vector normal to the surface, p(s) is the pressure, 
and C 5 (s) = z f (s)r n (s) - r'{s)z"(s) and C^s) = -z f {s)/r(s) are the principal radii of the curvature with 
primes (') denoting differentiation with respect to arclength. The unsteady Bernoulli equation may be used 
in (3.14) to relate the pressure to <£. This gives 

(ll) = ^«-^ 2 )-(C s + C,), (3.15) 

where the variables have been non-dimensionalized using the inertial scales 

length = R mass = pi? 3 time — (pR 3 /a eq ) 1 ^ 2 . (3.16) 

The time derivative in (3.15) is with respect to an observer moving with the normal velocity of the interface. 


dx 

~dl 

P 


3.2 Numerical procedure 

The double-layer potential boundary integral formulation uses equations (3.8), (3.9), (3.11), (3.12), (3.13) 
and (3.15) to follow the time-evolution of the drop’s shape. The equations are discretized by dividing the 
drop boundary into N equally-spaced nodes, with spacing As = L/(N - 1) in the interval 0 < s < L. 

Since the dynamics are driven solely by the curvature variations around the drop, it is important to 
calculate higher-order arclength derivatives accurately. For this purpose a least-squares spectral transform 
method [42] is used to represent all the parameterized surface variables. The method approximates each 
variable by a truncated sine and cosine series that is periodic in twice the total arclength L: 

„ x G 0 A f7Tks\ . f7Tks\ 

/(«) = y + ^ cos f — ] + ^Mmf — 1 

* k=l ' ' k=l V J 

2 m+ 1 

= X] c * S k{s). 

k= 1 

Here /(s) is any variable and Sk (s) is a sine or cosine shape function with 2 m + 1 = N c coefficients c * . 
For N discrete nodal values of /(s { ) and assuming N c < N, equation (3.18) represents an over-determined 


(3.17) 

(3.18) 


34 



system of equations for the coefficients. This may be solved in a least-squares sense [42]. The least-squares 
procedure results in a symmetric and positive-definite N c x N c system of equations for the coefficients c . *. 
Inverting the coefficient matrix is most efficiently accomplished using a Cholesky decomposition routine [36]. 
Once the function coefficients are determined, its derivative or integral may be exactly calculated. 

Aliasing is a potential problem in the least-squares spectral transform method. It can be avoided by 
adjusting N c so that the smallest wavelength represented by the shape functions, 4 L/(N C - 1), is larger than 
four times the largest of the unequally-spaced arclength increments As max that appear after each change in 
the drop’s shape. For the simulations shown in this work N c is chosen to be about half the total number of 
nodes N in order to accurately resolve the drop shape and eliminate aliasing. 

Since the drop is axisymmetric, all the surface variables are either even or odd functions in arclength. 
For instance, z and 4> are even functions while r and ip are odd functions. Higher-order derivatives of these 
functions are alternately odd or even functions. Neglecting the sine (or cosine) shape functions in (3.17) for 
even (or odd) variables decreases the size of the N c x N c system of equations for the least-squares spectral 
coefficients by about a factor of two. 

The inverse transform back to the nodal values of the variable involves a matrix multiplication by the 
discrete shape function matrix Sk(si). as defined in (3.18). To redistribute the nodal values of the variable at 
equal arclengths the discrete shape function matrix is evaluated at s* — (i — l)As,« = U ■ • ? N. Inspection 
of (3.17) shows that this matrix is independent of the total arclength L , which implies that it need only be 
calculated once for use in redistributing the variables. 

The regularized integral equations are discretized into matrix equations using a trapezoidal quadrature 
rule. For example, the discrete analog of the integral equation (3.8) takes the form 


N 

^ ^ Sj)^(Sj) — , for i 1, . . . , A , 

j= i 


(3-19) 


with 


K{si,Sj) 


f h{sj) ■ VG(5i,Sj)r(sj), for j / i, 

N 

1 - 5^(1 - 5 kl )w k n{s k ) ■ VG{si,s k )r{s k ) , for j = i. 
< k = 1 


(3.20) 


Here are the trapezoidal quadrature weights and Ski is the Kronecker delta symbol. The resulting 
matrix equation is solved using an LU decomposition routine. The use of regularized integral equations and 
trapezoidal quadrature weights allows for the complete discretization and solution of the integral equations 
without the need for explicit interpolations between the nodal values. 

The numerical procedure begins by initializing the shape of the drop x(s) and the scalar velocity potential 
<p(s) and calculating arclengths between the possibly unequally-spaced nodes. The arclengths between the 
nodes are first parameterized with the linear distance between the nodes. This approximation is then 
improved by: 1) solving for the least-squares spectral coefficients of r and z\ 2) calculating their derivatives 
r' and z ' with respect to the parameterization; and 3) integrating the arclength function [(r') 2 4- (z') 2 ] 1 / 2 to 
find the true arclength s at each node. 

The curvature, kernel, and tangential velocity at the interface is constructed from the derivative of the 
position vector and scalar velocity potential with respect to the true arclength. The normal velocity of the 
interface, however, must be found by solving the matrix equivalents of (3.8) and (3.9) with (3.11). Using 
the instantaneous velocity and curvature of the interface, the drop shape and scalar velocity potential are 
updated in time using a Runge-Kutta scheme. To prevent clustering, the nodes are redistributed to equal 
increments in arclength following each time step. 


3.2.1 Base flow: inviscid oscillations 

There are several conserved quantities that may be used to measure the accuracy of the numerical calculation. 
They include the volume, center of mass, and, for inviscid oscillations, the total energy, which is composed of 
the kinetic and potential energies. The quantities of volume and kinetic energy involve volume integrals, but 
these may be turned into surface integrals with the Divergence Theorem and the vector identity V • x = 3 
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in three dimensions. Since dS - r dsdip , these surface integrals may be turned into line integrals in the 
axisymmetric coordinate system. Thus, for the calculation of the volume of the drop, 

V = - f V • x dV = {zr l - rz*)rds . (3.21) 

3 J Vm 3 Jo 

Similarly, for the kinetic energy in the drop, since v n — n ■ ( V x^) = n-V</>: 

K.E. = f t:\v\ 2 dV = it f (j)v n rds . (3.22) 

Jv m 2 Jo 

The potential energy is the surface area of the drop 

P.E. = [ dS = 2tt [ rds. (3.23) 

Js, n 2o 

The initial conditions for the shape of the interface and the scalar velocity potential are taken to be 

|x(0)| = l+eP/(cos0) 

<f>{9) = 0 . (3.24) 

Figure 3.2 shows an example calculation for small-amplitude shape oscillations in the quadrupole mode. 
Calculations using a fourth-order Runge-Kutta time-stepping scheme did not differ significantly from the 
calculations using a second-order scheme. The second-order Runge-Kutta scheme was therefore used for 
all the calculations presented in this thesis. The calculated period of oscillations was found to be 2.22 in 
excellent agreement with the nondimensional theoretical value of 27r/V§ ~ 2.221. The conserved quantities 
of total energy and volume were found to deviate from their initial values by less than 10 7 % and 10 8 % 
respectively. 

3.3 Modified boundary conditions 

In this section, the base potential flow formulation for the drop in vacuum is modified to incorporate the 
effects of viscosity and surface rheology. Scaling the governing equations for this problem with the inertial 
scales, the nondimensional conservation of mass and momentum equations, and the transport equation for 
an insoluble surfactant take the forms: 




<1 

o 

= o, 

(3.25) 

Q w tot 

dt 

+ v tot 

. Vv tot 

= -Vp tot + a 2 V 2 v tot , 

(3.26) 

dY 

dt. 

+ v s - 

(v tot r) 

= 2— v s 2 r . 

Pe s 

(3.27) 


Here v tot and p tot are temporary notations for the nondimensional velocity and pressure fields. The surfactant 
transport equation is coupled to the continuity and momentum equations through the boundary conditions 
and the surface equation of state. Although a nonlinear surface equation of state could be incorporated into 
the numerical formulation, the linear surface equation of state 

<r(T) = l- e ;(r- 1) (3.28) 

is assumed here. This allows for direct comparisons to the corresponding results from the small-amplitude 
theory. The nondimensional components of the stress balance boundary conditions and the kinematic bound- 
ary condition take the respective forms 

p tot — 2a 2 (n ■ E tot • n) = (V, ■ n) - e^V, ■ n)(r - 1) 

+2 M ;(V S ■ E‘ ot ) • n + «;(V. • n)(V. • v tot ) 
a 2 (n ■ E tot • I s ) = -e*V,r + 2/i*(V s • E* ot ) • I s + k* s V s (V s ■ v tot ) 

tot \ 

— = n (n • v ‘) . 

dt 
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(3.29) 

(3.30) 

(3.31) 



Figure 3.2: Simulations of an initial perturbed in viscid drop in vacuum after the initial conditions (j> — 0 and 
|x(0)| = R[ 1 4- 0.01P 2 {cos6)] with N = 31, N c = 15, and At = 0.001. The scale factor R in the initial shape 
ensures that the nondimensional equilibrium volume V eq and the equilibrium energy E eq are 4 tt/ 3 and 47r, 
respectively. 





To begin to link the base potential flow formulation of the last section to the equations incorporating 
viscous and surfactant effects above, the total nondimensional velocity arid pressure are separated into two 
parts 

v tot = v + V (3.32) 

p tot - p + P. (3.33) 

The lower-case letters (v,p) represent the potential flow velocity and pressure fields while the upper-case let- 
ters (V, P) represent the vortical velocity and pressure fields. The potential flow fields satisfy the irrotational 
Euler equations 

0 (3.34) 

-Vp. (3.35) 

The vortical velocity and pressure fields satisfy the full Navier-Stokes equations (3.25) and (3.26), after 
subtracting (3.34) and (3.35): 

0 (3.36) 

-VP + a 2 V 2 v + a 2 V 2 V . (3.37) 

Inserting the decomposition (3.32) into the total surfactant transport equation (3.27) yields 
F)Y 1 

j- t + (v + V) • v s r + [v s ■ (v + v)]r = — v s 2 r . (3.38) 

The boundary conditions (3.29)-(3.31) may also be similarly decomposed, but for clarity this step is 
postponed for later sections. 

The boundary integral representation for the potential flow in the previous section drives this problem. 
The incorporation of the bulk viscous and surfactant effects comes in by replacing the boundary conditions 
for potential flow with the total stress balance and kinematic boundary conditions in the boundary integral 
formulation and keeping only the leading order terms that contribute bulk viscous and surfactant effects. 
The unknown vortical terms in the boundary conditions are found by satisfying the leading-order equations 
for the vortical fields. 

The vortical fields are assumed to be nonzero only in the thin Stokes boundary layers that form to satisfy 
the tangential stress balance boundary condition. This leads to considerable simplifications in subsequent 
sections. 

The idea is to use the order of magnitude estimates for the velocity and pressure fields from the small- 
amplitude theoretical analysis and keep all the terms in the modified boundary conditions (written in local 
coordinates) up to 0(a 2 ). These additional terms can be calculated to leading order from the Navier-Stokes 
equations (3.36) and (3.37). Doing so requires the assumption that the boundary layer fields decay away 
from the interface. 

This section begins by describing a method to include bulk viscosity in what would otherwise be a 
potential flow calculation for large Reynolds number shape oscillations of a drop in vacuum. The succeeding 
three cases include increasingly greater effects from the presence of an insoluble surfactant at the interface. 

3.3.1 Case 1: “negligible” surfactant effects 

This section considers the effects of bulk viscosity on the shape oscillations of a drop in vacuum when the 
nondimensional surface properties e*, p*, and k* are “negligible”, or 0(a 3 ). These effects are to be included 
to leading-order in the numerics by following the procedure just described. 

Order-of-magnitude analysis 

By using the relations in Appendix B, the nondimensional boundary conditions at the interface may be 
written in local coordinates as follows 

” +j ’- 2 “ i (^ + lf)= c ' +c - ,339) 


dv 


V -V = 

+ v ■ VV + V • Vv + V • VV = 


dv 

at 


V -V = 

+ v • Vv = 
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(3.40) 


a 


2 



fdV. dVn 

\ dn ds 



= 0, 


I - + (3,1) 

f t - !<* + V.). (3,2) 

Here r and 2 are the components of the interface shape. In (3.40), the tangential component of the stress 
balance, the symmetric part of the ns component of the rate of strain tensor corresponding to the potential 
flow part has been used. 

As before, it is convenient to introduce a new notation in the order of magnitude analysis. The small- 
amplitude analysis in section secxasel shows that the potential and vortical flow fields may be written 
as: 


v = n(u n0 + a 2 u n2 ) + s (u s0 + ol 2 u s2 ) + 0(a 3 ) (3.43) 

p = Po + ol 2 P 2 + 0(a 3 ) (3.44) 

V = na 2 t/ 7l2 + s (aUsi + ol 2 U s2 ) + 0(o?) (3.45) 

P = a 2 P 2 + 0(a 3 ). (3.46) 


As a reminder, the lower case fields (v,p) represent potential flow and the upper-case fields (V,P) represent 
vortical flow. Here the 0(a 2 ) vortical pressure P 2 has been introduced to account for the nonlinear pressure 
not appearing in the linear analysis. The potential and vortical fields depend on both n and s. The potential 
flow fields vary slowly in n and s over a nondimensional length scale of 0(1). The vortical fields, on the other 
hand, vary rapidly in n, and decay exponentially away from the interface over a nondimensional length scale 
of O(a). The normal derivatives of the vortical fields are denoted by a scaled normal coordinate n aN . 
For instance, the normal derivative of the 0(1) tangential component of the vortical velocity 


dU 8 o = 1 dUso 
dn a dN ’ 


(3.47) 


where dU s0 /dN is an 0(1) quantity. The order-of-magnitude estimates for the surfactant concentration T 
depend only on s and vary slowly over a nondimensional length scale 0(1). At this stage the only calculable 
quantities are the 0(1) quantities representing the base potential flow. 

Inserting the scales (3.43)-(3.46) into these boundary conditions and expanding yields 


Po + a 2 [ P 2 T P 2 - 2 nu ) — C s + + 0(q 3 ) 


d N 


2 ( dUnO _ + dU °' 


ds 


dN 


= 0 + 0(a 3 ) 


(3.48) 

(3.49) 


J = ~[u n0 +a 2 (u n2 +U n2 )}+O(a 3 ) (3.50) 

at os 

ft y dv 

^ = ^- s [u n0 + a 2 (u n2 +U n2 )] + O(a 3 ). (3.51) 

Neglecting the 0(ot 2 ) in these equations reduces them to those for the base potential flow. In order to 
incorporate the leading order effects of viscosity into the boundary integral formulation, the quantities at 
0(a 2 ) must be determined. Of these, the 0(a 2 ) corrections to the potential flows, u n2 and /; 2 , may be 
absorbed into the base flow calculation. Another term which may be treated at this stage is du no /dn, which 
may be calculated with the help of V • v = 0 written in local coordinates and evaluated at the interface 

dO nO du s o 1 dv (p i p \ n /q 

-air + -5“ 4- - -^ v s0 + (Cs + Cy)v n 0 = 0 . (3.52) 

dN ds r os 
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The three remaining vortical terms, U n 2 , P' 2 , and dU s \/dN , must be found from equations (3.36) and (3.37) 
for the vortical fields V and P. These fields need only satisfy the leading order equations since their higher 
order corrections are neglected in the modified boundary conditions above. 

From the assumed scaling, the continuity equation can be expanded to 0(a 2 ) 


V ■ V = a 


du n2 dUs 1 

1 ^ r 


dN 


ds 


;&«) - 


0(a 2 ). 


(3.53) 


The time derivative of the boundary layer velocity is 


av 

dt 


ds Tr . OU 


sl 


+ 0(a 2 ) . 


(3.54) 


Equation (3.54) involves the time derivative of the unit tangent vector. Since the unit tangent vector is the 
arclength derivative of the position vector, it can be related to the arclength derivative of the normal velocity 
through the kinematic boundary condition 


ds d f 9x 

dt dt V 9s 


d_ fdx 
ds \ dt 


u n o) + 0{a ) . 


It follows that 


dV . du n o 


) = l ( " 

^ T CsUnoUsl'j + C^( a2 ) • 


The nonlinear terms in (3.37) may be also expanded to give 

dU n 2 


V • VV = na ( u n Q 
+s 


^nO 


dN 

dU sl 


) 


V ■ Vv — n a 


tSft 


dN 
du n Q 


^s^sO P s\ 

9Us2 

dN 


+ Ot u n0 


+ w s0 * 


dU s 


ds 


+ 0(a 2 ) 


ds 
/ du s o 
V ds 


Cs^sO ) bsl 


+ C s u n o ) U s i + 0{a 2 ) 


V ■ VV = 0{a 2 ) . 

Similarly, for the pressure and viscous terms 


d P 

VP = na^ + 0(a 2 ) 


a 2 V 2 v = 0(a 2 


(3.55) 


(3.56) 


(3.57) 

(3.58) 

(3.59) 

(3.60) 

(3.61) 


a 2 V 2 V = s Q ^ + 0(a 2 ). 

dN z 


Combining (3.54) — (3.62) in the Navier-Stokes equations gives, to leading order: 

+ 




and 


dU s i 

OLU 7 iO ® 


f 9lL n Q 
V ds 


■c 


du sl , 9U.1 , du s2 , dU sX , _ 

a— — — b u no~rncr + au n0 — r + <xu s o— I- a 


\ 1 dr 

•'S^'sO ) o ^riO 

) T ds 
( du s0 


dP 

Us\ = a ^K Q ~) 


at 


aw 


aw 


as 


^ ^ + 2C s u n0 ^ 17*1 — a Qpj 2 + ) 


(3.62) 

(3.63) 

(3.64) 

(3.65) 
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The continuity equation (3.53) may be identically satisfied by defining the vortical part of the boundary 
layer velocity as the curl of a vector velocity potential 'J' = In general V takes the form 

(3.66) 


V = Vx$ 


, d h s dh v _ 

hs^r + T-r* 

(sS ds 


+ s 


— "'ll 

on 


fan dh^ 

On 


* 


(3.67) 


Here, h n , h s , and are the metrical coefficients for the local coordinate system. To obtain order of 
magnitude estimates, the vector velocity potential is written $ = a 2 9 2 + O(a^), with ^2 assumed to have 
large normal gradients. The velocity magnitudes may be thus restated in terms of $2 so that 


U n 2 = h. 


Us 1 

U S 2 


— fan 


d^2 hs dh, 

ds 

d^2 


h w ds 


^ 2 


dN 


fan dhip 


(3.68) 

(3.69) 

(3.70) 


With the continuity equation identically satisfied, these expressions may be substituted into equations 
(3.63) — (3.65) for the vortical fields. Notice that with these substitutions every term in the equations con- 
tains the vector velocity potential, which decays exponentially away from the interface with a decay factor 
1 /a, and at least one normal derivative. These observations, and the knowledge that all the other quantities 
in the equations vary slowly over the thin Stokes boundary, allow for the simplification of the equations 
through an integration in the normal direction. In a sample term containing an explicit normal derivative, 
this process proceeds by first expanding the slowly varying terms about their values at the interface, 


mrno^ = a«.„«U 0 ^ + O(a 2 ) (3-71) 

and, by integrating over a nondimensional boundary layer thickness a, where the nondirnensional variable n 
is rescaled with a, inside the drop 

dU sl 


f 

J — tx 


dU *l ^AJ 

aUn0 ~dN adN 


* w no| n =o / 

J — ex 


dN 


dN + 0(a s ) 


Unc\ n=0 (U sl \ n = 0 -U sl \ n= _J+O(a 3 ) 

d$ 2 ' 


WnO 


dN 


+ 0{a i ) 


(3.72) 


n=0 


where the last expression on the right-hand side has been evaluated at the interface n — 0, where the metric 
coefficients are known. 

Other terms, may be similarly approximated: 


du s o 

a—z—Usi = 


ds 


du s0 , d ^ 2 

Ot fan' 


= —a 


— —a 


ds 
du s o 


ds 

duso 


dN 

n = 0 ^ 


ds 


n=0 


, d* 2 \ 

hn dN ) 


+ 0(a 2 ) 

dh n 


a — ^2 

on 


+ O(or 


(3.73) 


The second term in the square brackets may be neglected. This expression may now be integrated over the 
nondimensional boundary layer thickness to provide 




du s o 
ds 


U s \ adN — —a" 


du 


sO 


d 


— —a 


— -a 


ds 

dugQ 


f ^rj{h n ^ 2 )dn + 0(a 3 ) 

n=0 J -oo 


ln=0 


/ dUsQ 
\ ds 




[(^2)l n=0 - (hn* 2 )\ n= - 1 ]+ 0 (a 3 ) 

+ 0(a 3 ) 


(3.74) 


n=0 
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(3.75) 


Proceeding similarly for each of the terms in the vortical equations leads to, 


a 2 P2 = —OL 2 u n o~q^ + 


du n 0 \ 1 9/* 

21 c - 5 " s0 j ras u,, ° 


$2 + 0{a z 


, 5^2 < 9^2 

° ~df + aUn0 dN 


(3.76) 


All the terms up to 0(a 2 ) are retained and evaluated at the interface. The remaining normal derivatives in 
these equations can be evaluated as follows. The normal derivative on the left side of (3.76) can be absorbed 
into the time-derivative of the vector velocity potential with respect to an observer moving at the velocity 
u n0 . The normal derivatives on the right sides of the (3.76) may be found from the tangential stress balance 
boundary condition. 


Modified equations and boundary conditions 

Dropping the order of magnitude estimates and reverting to the original notation in terms of v, V, and $ 
yields the following modified boundary conditions for the base potential flow 


where 


~Dt j 
DV \ 
~Dt)n = 

dr 

dt 

dz 

dt 


(v 2 n - vj) - (C, + C v ) - 2 « 2 ^ + v n V„ + P 


-V s 


CM 

ds 


dv 

— 1- (2 C s T C^)v n 

os 


+ a : 


dri 2 


dz 

~d~s {Vn + K) 
%(v n + V n ) 


(3.77) 

(3.78) 

(3.79) 

(3.80) 
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2 dv n 
dn 

V n 

P 

,d 2 fr 

dn 2 


—a 

d<n 


dv s 1 dr 

1 — 7T Us s + C^)v n 

os r os 


1 dr 

_1 ip 

ds r ds 
d ^ 

~ v "-g; + 


= -2 a 2 


dv n 

ds 


( dv n \ 19r ] 

2 -5 C s v s - ~jr v n * 

\ os J r os 

-C s v^j 


(3.81) 

(3.82) 

(3.83) 

(3.84) 


These equations are evaluated at the interface and retain all the leading order viscous terms at (9(a 2 ), but 
neglect the higher order contributions. Note that in (3.77), Bernoulli’s equation has been used to relate 
the pressure to the time derivative of the scalar velocity potential and a V n (d(j>/dn) = v n V n , which is an 
0(a 2 ) term, has been added to both sides to ensure that the time derivative is consistent with the kinematic 
boundary condition. A similar term V n (d^/dn) has not been added to (3.78) because it is an 0(a 3 ) term 
and may be neglected in the above equations. 

The numerical procedure begins by initializing the r(s), z(s), and 'P ( .s ) . After calculating the 

arclengths, the kernels of the boundary integral formulation are formed and used to calculate the base 
potential flow velocity components. The potential flow velocity components are used to evaluate the right- 
hand sides of the evolution equations for r, z, 0, and tp, which are updated with the Rung-Kutta time 
stepping scheme. 

Figure 3.3 shows an example calculation for small-amplitude quadrupole shape oscillations of a weakly 
viscous drop. The additional evolution equation and modified boundary conditions were seen to produce an 
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Figure 3.3: Simulations of an weakly viscous drop in vacuum for a 2 = 0 (dashed line) and a 2 = 0.005 
(solid line) after the initial conditions 4> = 0 and |x(0)| = i?[l + 0 . 01 P 2 (cos ^)] with N = 31, N c = 15, and 
At = 0.001. The scale factor R in the initial shape ensures that the nondimensional equilibrium volume V eq 
and the equilibrium energy E eq are 4n/3 and 47T, respectively. 


instability in the simulations. The instability was traced to the variable V, which tended to develop small 
wavelength disturbances that grew in time. Since V and its first derivatives in arclength are used to calculate 
Vm these high wavenumber oscillations eventually contaminated all the variables and the calculations had 
to be stopped. The calculations were found to be stable in time if the V n were eliminated from the evolution 
equations for the shape. The damping constant, which results from the potential flow term —2 a 2 dv n /dn, 
did not correspond to that predicted from theory and, interestingly, was off by a factor of 2 for all axi sym- 
metric mode numbers. The instability arising in the vortical velocity potential V has been well documented 
in previous work (see, for example [26]). The customary method of correction introduces a “smoothing” 
parameter D into the problem, where the nodal values of ’F were replaced with 

V { = Vi - D(Vi-2 - 4^-1 + 6Vi - 4 Vi-i + V i+2 ) (3.85) 

after each time step. This process amounts to placing the negative of a fourth-order derivative with respect to 
arclength on the right-hand side of the evolution equation for and acts to smooth out the short wavelength 
disturbances appearing in the variable. 

Preliminary results proved that, with an appropriate choice of the parameter D , the process of smoothing 
eliminated the instabilities. It was noticed, however, that with a slightly different smoothing algorithm 

Vi = Vi + D{Vi - 2 + 16V i-i - 30 Vi + 16^-1 - * i+2 ) , (3.86) 
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Figure 3.4: Simulations of a weakly viscous drop in vacuum for a 2 = 0 (dashed line) and a 2 = 0.005 (solid 
line) after the initial conditions 0 = 0 and |x(0)| = R[ 1 4- 0.01P2(cos #)] with N = 31, N c = 15, and At = 
0.001. The scale factor R in the initial shape ensures that the nondimensional equilibrium energy E eq is 47 t. 
Also shown is the theoretical prediction for the damped shape in this limit z iop (t) = 1 + [zt 0l) (t = 0) - l]e“ 5a 1 
(dotted line). 


the size of the smoothing parameter D could be decreased by two orders of magnitude. This different 
algorithm amounts to placing a second-order derivative with respect to arclength on the right-hand side 
of the evolution equation for V P. If the smoothing parameter was set too large, the damping constant was 
less than the theoretical value. A smoothing parameter of D — 10 -4 was used to produce the weakly 
viscous simulations in Figure 3.4. With smoothing, the damping rate for the shape oscillations was in good 
agreement with the theory. As in the base potential flow calculation, the volume was found to deviate by 
less than 10~ 8 % from its initial value. 

3.3.2 Case 2: “small” surfactant effects 

This section considers the effects of bulk and surface viscosity and surfactant transport on the shape oscilla- 
tions of a drop in vacuum when the nondimensional surface properties e*, ^*, and k* s are “small”, or 0(a 2 ). 
These effects are to be included to leading-order in the numerics by expanding the total decomposition of 
the boundary conditions in local coordinates to 0(a 2 ), using the order-of-magnitude estimates from the 
corresponding section in the theory, and neglecting those terms appearing at 0(a**) and higher. 

Order-of-magnitude analysis 

With the velocity and pressure decomposed into potential and vortical fields, the nondimensional forms for 
the normal stress balance, tangential stress balance, and kinematic boundary conditions in local coordinate 
are, respectively, 

p + P-2a 2 (^ + ^f) =(C s +C VJ )-e:(C s +C^)(r-l) 

+A Ci.Cs ~ C v )(m~ + M~) + k*{C s + C v ){m + + M + ) , (3.87) 
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(^n "b V n ) , 


where 


M ± = 


^±~«. + (C.±C>„) 

as r as / 

^±l|:v. + «:.±c,)v„ 

as r as 


71—0 


71 = 0 


(3.88) 

(3.89) 

(3.90) 

(3.91) 

(3.92) 


The order-of-magnitude estimates for the potential and vortical flow fields and the surfactant concentra- 
tion from the corresponding section in the theory yield expansions of the form 


V 

= n (u n o + a 2 u„ 2 ) + s (u s0 + a 2 u s2 ) + 0(a 3 ) 

(3.93) 

p 

= po + a 2 p 2 + 0(a 3 ) 

(3.94) 

V 

= na 2 U n i + s (aU s i + a 2 U s2 ) + 0{o?) 

(3.95) 

P 

= a 2 P 2 + 0(a 3 ) 

(3.96) 

r 

= l + r 0 + ar\ +a 2 r 2 + 0(a 3 ), 

(3.97) 


where, as before, the lower case fields (v,p) refer to those resembling potential flow and the upper-case 
fields (V, P) refer to those resembling vortical flow. Here the 0(a) vortical pressure Pi has been introduced 
to account for the nonlinear pressure not appearing in the linear analysis. The potential and vortical 
fields depend on both n and s, but unlike the potential flow fields, which vary slowly in n and s over a 
nondimensional length scale of 0(1), the vortical fields vary rapidly in n, decaying exponentially away from 
the interface over a nondimensional length scale of O(a) and slowly in s. The normal derivatives of the 
vortical fields are denoted by a scaled normal coordinate n — > aN. For instance, the normal derivative of 
the 0(1) tangential component of the vortical velocity 


dU s0 

dn 


ldU s0 

a dN 


(3.98) 


where dU s o/dN is an 0(1) quantity. The order-of-magnitude estimates for the surfactant concentration T 
depend only on s and vary slowly over a nondimensional length scale 0(1). 

Substituting the order-of-magnitude scales (3.93)-(3.97) into the boundary conditions (3.93)-(3.97) and 
expanding to 0(a 2 ), noting that a normal derivative of a vortical field introduces a factor of 1/a, yields 


Po T at ( p *2 + P ‘2 ~ 2- 


du 


TlO 


dn 


{c, + c v )-e;(c. + c v ){ r 0 -i) 


+/i* (C s - C^)m 0 + k* s (C s + C v )m^ + ©(a 2 ) , 


(3.99) 
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(3.100) 
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and 


% = -^[Un 0 +a 2 (Un 2 +U n2 )} + O(a 3 ) (3.101) 

dt ds 

J t = ^ s [u n0 + a 2 (u n2 + U n2 )} + O(a 3 ), (3.102) 


where 




du s0 

ds 


1 dr 
r ds 


n s o H“ (^s di C^)u n o 


n=0 


(3.103) 


If the 0{a 2 ) terms are neglected in these equations, they simplify to those for the base potential flow, 
provided the unsteady Bernoulli equation is used to relate the pressure to the 0{ 1) scalar velocity potential. 
In order to incorporate the leading-order effects of bulk and surface viscosity and surfactant transport into 
the boundary integral formulation, the quantities at 0(a 2 ) must be determined. Of these, u n 2 and p 2 may 
again be absorbed into the base flow calculation and du no ldn calculated from Laplace’s equation written 
in local coordinates and evaluated at the interface. Note that apart from the e* terms all the quantities 
associated with the surface properties are known from the potential flow and may be readily evaluated. 

The unknown quantities still to be determined include the P 2 , dU s i/dN , and the tw r o e* terms. Of 
these P 2 and dU s 1 /dN may be found by introducing a vector velocity potential for the vortical velocity and 
applying an or der-of- magnitude analysis to the equations for the vortical fields to determine an evolution 
equation for the vector velocity potential and an expression for the vortical pressure in terms of this vector 
velocity potential both evaluated at the interface, as in the previous section. 

For the surfactant concentration, an evolution equation may be found from an order of magnitude analysis 
of the decomposed surfactant transport equation. Inserting the scales into this equation gives, to leading- 
order, 
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or combining terms and rearranging, 
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(3.104) 

(3.105) 

(3.106) 

(3.107) 


(3.108) 


Here the time derivative may be rewritten as 


ar 0 

dt 



n 


since the insoluble surfactant monolayer has no gradients in the normal direction. 


(3.109) 


Modified equations and boundary conditions 

Dropping the order of magnitude estimates and reverting to the original notation in terms of v, V, and 
T yields the following modified boundary conditions for the base potential flow 
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These equations are evaluated at the interface and retain all the leading-order bulk and surface viscous and 
surfactant transport terms at O (a 2 ). 

The numerical procedure begins by initializing the r(s), z(s), 0(s), \&(s), and T(s). After calculating 
the arclengths, the kernels of the boundary integral formulation are formed and used to calculate the base 
potential flow velocity components. The potential flow velocity components are used to evaluate the right- 
hand sides of the evolution equations for r, 2 , </>, 4', and T, which are updated with the Rung-Kutta time 
stepping scheme. 

The simulations including weak bulk viscous and surfactant effects were even more susceptible to the 
kind of instabilities in the variable ^ seen in the previous case including weak bulk viscous effects. No 
value of the smoothing parameter could be found for stable and physically realistic simulations including the 
surface viscosities and Gibbs elasticity that were of the same order as the inverse of the Reynolds number. 
For instance, with the Gibbs elasticity set to zero and both the surface viscosities nonzero, the smoothing 
parameter had to be set so large that the resulting oscillations damped out at a rate smaller than the case 
without surface viscosities. 

There were tw T o limits in which stable and physically realistic simulations could be obtained. These limits 
set the Gibbs elasticity and one of the surface viscosities to zero. In both example calculations the smoothing 
parameter was increased to D = 7.5 x 10“ 4 . Figure 3.5 shows the case when the surface shear viscosity is 
set to zero, but the surface dilatational viscosity is the same order as the inverse of the Reynolds number. 
The damping constant for small oscillations was found to agree with the theory in that limit. Figure 3.6 
shows the case when the surface dilatational viscosity is set to zero, but the surface shear viscosity is the 
same order as the inverse of the Reynolds number. The damping constant for small oscillations was found 
to agree with the theory in this limit. For the quadrupole mode, the shear viscosity was seen to have a much 
larger effect on the damping of the shape oscillations. 


3.3.3 Case 3: “medium” surfactant effects 

This section considers the effects of bulk and surface viscosity and surfactant transport on the shape oscil- 
lations of a drop in vacuum when the nondimensional surface properties e*, //*, and k* s are “medium”, or 
0(a). These effects are to be included to leading-order in the numerics by expanding the total decomposition 
of the boundary conditions in local coordinates to 0(a), using the order-of-magnitude estimates from the 
corresponding section in the theory, and neglecting those terms appearing at 0(a 2 ) and higher. 
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Time 


Figure 3.5: Simulations of an initially perturbed weakly viscous drop in vacuum with small surfactant effects 
for (a 2 ,e*,/is,K*) = 0 (dashed line) and (a 2 ,/s*) = 0.005 with (e*,K*,jx*) = 0 (solid line). The initial 
conditions are (0, 'F) = 0 and |x(#)| = P[1 + 0.01p2(cos#)] with N — 31, N c = 15, and At = 0.001. The 


scale factor R in the initial shape ensures that the nondimensional equilibrium energy E eq is 47 t. Also shown 
is the theoretical prediction for the damped shape in this limit z top (t) = 1 + [zt op (t = 0) — l]e _ ^ 5a 


(dotted line). 
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Figure 3.6: Simulations of an initially perturbed weakly viscous drop in vacuum with small surfactant effects 
for (a 2 ,e*,/i*,/c*,Pe 5 ) = 0 (dashed line) and (a 2 ,/z*) = 0.005 with (e*, /c* T Pe s ) = 0 (solid line). The initial 
conditions are (</>, 4') = 0 and |x(0)| = i?[l + 0.01P2(c°s^)] with N = 31, N c = 15, and At — 0.001. The 
scale factor R in the initial shape ensures that the nondimensional equilibrium energy E eq is 4n. Also shown 
is the theoretical prediction for the damped shape in this limit z t op(0 = 1 + [ztop(t = 0) — l]e“^ 5a 
(dotted line). 
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Order-of-magnitude analysis 

This section applies an order of magnitude analysis to the general expressions for the normal and tangential 
stress balance and kinematic boundary conditions, the equations for the vortical fields and the surfactant 
transport equation, written in local coordinates. 

The order-of-magnitude estimates for the potential and vortical flow fields and the surfactant concentra- 
tion, obtainted from the corresponding section in the theory, are given by 


V 

= h(u n o + att nl ) -f s (u s q + au s i) -t- 0(a 2 ) 

(3.120) 

p 

= po + api + 0(a 2 ) 

(3.121) 

V 

= haU n \ + s (U$o + ol 2 U s \ ) + 0{a 2 ) 

(3.122) 

P 

— aP\ -f (9(a 2 ) 

(3.123) 

r 

— 1 + Tq + aT i + O (a 2 ) , 

(3.124) 


where, as before, the lower case fields (v,p) refer to those resembling potential flow and the upper-case 
fields (V, P) refer to those resembling vortical flow. Here the 0(a) vortical pressure Pi has been introduced 
to account for the nonlinear pressure not appearing in the linear analysis. The potential and vortical 
fields depend on both n and s, but unlike the potential flow fields, which vary slowly in n and s over a 
nondimensional length scale of 0(1), the vortical fields vary rapidly in n, decaying exponentially away from 
the interface over a nondimensional length scale of O(a) and slowly in s. The normal derivatives of the 
vortical fields are denoted by a scaled normal coordinate n — > aN. For instance, the normal derivative of 
the 0(1) tangential component of the vortical velocity 

dUso _ IdUso (3.125) 

dn a dN ’ 

where dU s0 /dN is an 0(1) quantity. The order-of-magnitude estimates for the surfactant concentration T 
depend only on s and vary slowdy over a nondimensional length scale 0(1). 

Substituting the order-of-magnitude scales (3.120)-(3.124) into the boundary conditions (3.93) (3.97) 
and expanding to 0(a 2 ), noting that a normal derivative of a vortical field introduces a factor of 1/a, yields 

Po + a ( pi 4* Pi ) = (C s + Cp) — e*(C s + C v? )To 

+/i;(C, - Cf)(m o + Mo) + <(C S + 0,)(m+ + M+) + 0(tx 2 ) , (3.126) 


and 


where 


si 


da 

* dN 

+/f 


= — e 


dT 0 

ds 


d , 2 dr _ 

— (m 0 + M 0 ) + + M o 


+ K sQ^( m o + M o) + 0( a2 ) i 


rjj' (Q y 

= — [u n o + a(u n i + U n \ )] + 0(a 2 ) 

-J7 = jr-[ u nO + a ( u nl + Uni)] + 0(ot 2 ) , 

at os 


™ o 


du s0 1 dr 

0 i r\ ^sO (C s i C^)u n o 
ds r ds 

dU s o , 1 dr 

r\ -k r\ b SO 

ds r os , , n=0 


ln=0 


(3.127) 

(3.128) 

(3.129) 

(3.130) 

(3.131) 


If the 0(a) terms are neglected in these equations, they simplify to those for the base potential flow, provided 
the unsteady Bernoulli equation is used to relate the pressure to the 0(1) scalar velocity potential. In order to 
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incorporate the leading-order effects of bulk and surface viscosity and surfactant transport into the boundary 
integral formulation, the quantities at 0(a 2 ) must be determined. Of these the p\ may again be absorbed 
into the base flow calculation. 

The unknown quantities still to be determined include the Pi, To, dV s o/dN , and the two e* terms. Of 
these Pj and dV s0 /dn may be found by introducing a vector velocity potential for the vortical velocity and 
applying an order-of-magnitude analysis to the equations for the vortical fields to determine an evolution 
equation for the vector velocity potential and an expression for the vortical pressure in terms of this vector 
velocity potential both evaluated at the interface, as in the previous section. The U s o and its arclength 
derivatives, however, cannot be determined using the techniques of previous sections. This is because the 
tangential component of the vortical velocity is defined in terms of the normal derivative of the vector velocity 
potential, which cannot be determined using any variable known only at the interface. Without the ability 
to calculate the tangential component of the vortical velocity explicitly, the techniques used in the previous 
cases may not be used to calculate the numerics for case 3. 

Boulton-Stone [12] has developed alternate techniques to solve for the tangential velocity of the interface 
explicitly, by generating a thin grid near the interface which remains locally orthogonal in the local coordinate 
system, but this analysis is not attempted in this work. 


3.3.4 Case 4: “large” surfactant effects 

This section considers the effects of surface viscosity and surfactant transport on the shape oscillations of 
a drop in vacuum when the nondimensional surface properties e*, /i*, and k* are “large”, or greater than 
O(a). The formulation for this case neglects all effects at 0{a) and smaller. 

For this case the bulk viscous dissipation and the effects of the vortical fields may be neglected to leading 
order. The flow field is calculated from potential flow equations and the shape of the interface is updated 
with the normal velocity from potential flow. Since the bulk viscous forces are neglected the tangential 
component of the interface velocity, however, is different from the tangential velocity from potential flow. 

The governing equations are the regularized integral equations representing Laplace’s equation for the 
potential flow in the drop and the surfactant transport equation. These equations are subject to the normal 
and tangential components of the leading-order stress balance boundary conditions 

P = (C S + c v ) - e ;(c s + c„)( r - 1 ) + p s (c s - c v )m~ + k‘ s (c s + cjm+ (3.132) 


°-~ e; ^ + u: 


dM- 

ds 


2 dr 

H — 

r os 




dM+ 

ds 


and the leading-order kinematic boundary condition 


+ C*(u n i + U n \ )] + G(a 2 ) 

— ^[w n 0 + Ct(u n l + Uni )] + 0(a 2 ) , 


(3.133) 


(3.134) 

(3.135) 


In the stress balance boundary conditions, and in the surfactant transport equation, 


M ± = 



(3.136) 


and is composed of the normal and tangential velocity components of the interface, where the normal interface 
velocity is the normal component of the potential flow velocity evaluated at the interface. 

The tangential stress balance equation may be written in the form of a second-order differential equation 
for the unknown tangential interface velocity 


;) 2 v ft V 

a{s)~~{s) + 6(a)^(s) + c(«)V,(«) = d{s) (3.137) 

where a, 5, c, and d are known functions from the potential flow and the shape. Since the tangential interface 
velocity is necessarily and odd function in arclength for this axisymmetric problem, it may be expanded in 
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a truncated sine series 


n c 

v*(s) = Y J Ck sin 

k=l 



(3.138) 


where N c is the number of coefficients limited by the Nyquist frequency. For discrete nodal values Si,i — 
1 . . . N in arclength the tangential stress balance boundary condition may be written 


where 



(3.139) 


(3.140) 


If N > iV c this equation represents an over- determined system of equations for the coefficients of the tan- 
gential interface velocity, and may be solved in a least-squares sense by multiplying each side of the equation 
by the transpose of A **. The tangential stress balance boundary condition becomes a matrix equation for 
the coefficients of the tangential interface velocity. 

The modified boundary conditions for this case are 


D<f>\ 

Dt) n 

~ — (C s + ^ip) + e;(C s 4- C*)(r 


-n*(c t - c v )M~ - k.;(c, + c v )m + 

d r\ 
Dt) n ' 

1 d 2 T dr + 

Pe s ds 2 Vs ds 

dr 

dz 

dt 

— h - 
as 

dz 

dr 

dt 

Js Vn ‘ 


(3.141) 

(3.142) 

(3.143) 

(3.144) 


where the discrete form of the tangential stress balance boundary condition (3.140) is used to determine the 
tangential velocity of the interface. 

The numerical procedure begins by initializing the r(s), z(s), 0(s), and T(s). After calculating the ar- 
clengths, the kernels of the boundary integral formulation are formed and used to calculate the base potential 
flow velocity components. Using the known interface shape and base potential flow velocity components and 
the coefficients of the tangential velocity of the interface are found from. Finally, the right-hand sides of the 
above evolution equations for r, z, 0, \I>, and T are updated with the Rung-Kutta time stepping scheme. 

Without the need to calculate the vector velocity potential for the vortical velocity contribution, the 
simulations for this case suffered none of the numerical instabilities seen in the previous cases. The procedure 
using the tangential stress balance to solve for the tangential velocity of the interface worked well and allowed 
for the combinations of surface viscous and Marangoni effects to be analyzed. 

Figure 3.7 shows a calculation in the limit when the Gibbs elasticity and surface Peclet number are zero 
with nonzero surface viscosities. Other simulations in this limit with 0(1) surface viscosities showed an 
over-damped behavior. 
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Figure 3.7: Simulations of an initially perturbed inviscid drop in vacuum with large surfactant effects for 
(e* , k * , fj ,* , Pe s ) = 0 (dashed line) and (*;*,//*) = 0.05 with (e*,Pe 5 ) = 0 (solid line). The initial conditions 
are (<j>, ’F) = 0 and |x(0)| = R[ 1 4- 0.01 /*2 (cos 0)] with N = 31, N c = 15, and At = 0.001. The scale factor 
R in the initial shape ensures that the nondimensional equilibrium energy E cq is 47r. Also shown is the 
theoretical prediction for the damped shape in this limit ztop(^) = 14* [zt op (t = 0) — l]e”l 16,VjMj ^ 3h '* +2/lj ^ 
(dotted line). 
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Chapter 4 

Summary and conclusions 


In this research, the dynamics of a drop in a fluid medium with a surfactant contaminated interface are 
analyzed. The fluids are modeled as Newtonian, and the insoluble surfactant monolayer is modeled as a 
Boussinesq-Scriven Newtonian surface fluid. Using a linear relationship between the local surface tension 
and surfactant concentration, an exact mechanical energy equation was derived for the system. The resulting 
nondimension al total mechanical energy equation is 

— {K.E. + P.E.} = — c* 2 {Bulk Diss.} 

dt 

-/i* {Surf. Diss.} 

+e*{Marangoni} . (4.1) 


The terms in (4.1) are clearly identifiable as volume and surface integrals over each phase. In particular, 
the physical interpretation of the Marangoni term was shown to be a combination of energy storage and 
dissipation terms whose dominant contributions are controlled by the surface Peclet number. 

A novel analysis of the total mechanical energy equation was performed using an averaging method. This 
method, when complemented with velocity field approximations found using matched asymptotic expansions, 
allows for the construction of a damped harmonic oscillator equation that approximately describes the 
dynamics of the system for small-amplitude shape oscillations. For the cases with surfactants this damped 
harmonic oscillator equation was coupled to a second o.d.e. that describes the surfactant transport. 

For inviscid potential flow and a clean interface, the dynamics of the system are approximately described 
by the simple harmonic oscillator equation: 


Ci£ + 


til- 1)(* +!)(* + 2) 

[(«+!)+ tp/p] 1 


(4.2) 


The nondimensional frequencies predicted by (4.2) are in agreement with the classical results of Lamb [23]. 

Adding viscosity to the system with a clean interface introduces damping. As shown in section 2.4.2 
(Case 1), the shape oscillations of the viscous system are described by the equation 


(1 + a Afi) &£ + (aBfi + a 2 B^2) &£ -h D 2 0 a t ~~ 0 • 


(4.3) 


Equation (4.3) differs from (4.2) in that it contains added mass and damping terms due to the bulk viscosity 
in both phases. The frequencies and damping constants for this case are accurate to 0{a z ) and agree with 
the results of Marston [28]. 

Even small surfactant effects lead to qualitative changes in the system. As shown in section 2.4.3 (Case 
2), for a viscous drop in vacuum with weak 0(a 2 ) surfactant effects the dynamics of the system are described 
by a coupled set of equations: 

at + at [o?2{t - l){2t + 1) + K* a t(t - l) 2 + - l)(t + l)(£ + 2)] 

+ at[l{l -!)(< + 2)] = - 5< [e;^-l)], (4.4) 
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9t + ~ 9e = (l - \)at . (4.5) 

In this case, the size of the surface Peclet number Pe a controls the modification of the frequency and damping 
constant. 

The dynamics of the system becomes more complicated for slightly larger (“medium”) surfactant effects. 
Section 2.4.4 (Case 3) shows that when k* and /r* are 0(a), a damped harmonic oscillator equation can be 
found only when the Gibbs elasticity is zero. For the viscous drop in vacuum, that equation is 


+ 


+ 


at { 1 + 


(*+l) 


s/n^e(2e + 1) 


\Ceol 21 

y/2 {2£+l) 


(Ceoe 


-17t/4 


+ c.c.) | 


at ^ ot \/ flro ^ ^ \Ceo \ 2 


+H* S (£ - l)(f + 1)(^ + 2) [1 + (C e o + c.c.) + \C 0Q \ 2 ] 

+k*£ [(£ — l) 2 + {£ — 1)(^ + l)(Cflo + c.c.) + {£ + l)-|Cflo| 2 ] | 
a t [£{£- !)(£ + 2)] = 0 


(4.6) 


The added mass and damping terms in this equation are relatively complicated functions of the Reynolds 
number and surface viscosities. 

More straightforward results are available when the surfactant effects are 0(1) and completely dominate 
the bulk viscous effects in the drop/mediurn system (Case 4). For that case, the following coupled system 
of equations is obtained in section 2.4.5: 


4- 1) + £p/p] + 

+ 


. f — l)(l + l)(l + 2) \ 

af \[^+l)+K(^-l)^ + 2)]J 
atW-m+ l)(^ + 2)] 

( 2cX^-l)(f+l)(f + 2) | 

9e \{K' s e(e + i) + fii(e-w + 2)}t 


9t + 


f 1(1 + 1) , e;^+l) 

1 Pe s + [K*f(f+l) + p*(f- 1)^ + 2)] 


-at 


2^(£ -l)(l + 2) 

[«;/(< +1)+MI(« -!)(< + 2)] 


} 


(4.7) 


(4.8) 


Supplementing the small-amplitude theory just summarized, a numerical boundary integral equation 
formulation was developed for large-amplitude shape oscillations of a drop in vacuum. With an asymptotic 
analysis of the weakly singular integrands in the regularized integral equations, a discrete formulation for 
their solution was developed that only used the nodal values of the variables without special treatment near 
the singularities. To accurately calculate the higher-order arclength derivatives in the formulation, a least- 
squares spectral transform technique was used. The least-squares representation allowed for clean coding and 
gave a continuous, smooth representation of the dependent variables. Such a representation was essential 
in evaluating derivatives, integrals, and arclengths between unevenly-spaced nodes. It further provided an 
efficient way of solving for the tangential velocity of the interface when surface parameters were large. The 
results for the base potential flow calculations were validated with excellent agreement with theory for small- 
amplitude oscillations. With a second-order Runge-Kutta. time-stepping scheme, example calculations for 
small-amplitude oscillations in the quadrupole mode conserved the total energy and volume with errors of 
less than 10 -7 %. 

Using the order-of-magnitude estimates for the velocity components from the theory, the boundary in- 
tegral equation formulation for potential flow was extended to approximately include the effects of bulk 
viscosity, surface viscosity, and surfactant transport. These effects were approximated by including only the 
dominant terms (in a) in the general viscous boundary conditions at a surfactant-laden interface. For the 
case including bulk viscous effects with a clean interface, an instability was seen in the calculations. This 
well-known numerical instability [26, 11, 44] could be controlled with a five-point smoothing algorithm. The 
resulting calculations gave damping constants that were in excellent agreement with small-amplitude theory. 
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Section 3.3.2 examined bulk viscous effects and small surfactant effects in the numerics. In this case 
the numerical instability was much more difficult to control. Only the limit when either one of the surface 
viscosities was nonzero could be calculated with the smoothing algorithm. The damping constants calcu- 
lated for each limit were in agreement with small-amplitude theory. For small-amplitude oscillations in the 
quadrupole mode, the additional damping from surface shear viscosity alone was an order of magnitude 
larger than that from surface dilatational viscosity. 

Case 3 included bulk viscous and surfactant effects assuming the surface properties to be 0{a ). The 
order-of-magnitude analysis for this case revealed that the tangential component of the vortical velocity was 
needed explicitly in the modified boundary conditions. Since this quantity is not available in the boundary 
integral formulation presented here, this case could not be treated. 

When the surfactant effects were large (Case 4), the vortical velocity contribution could be neglected 
entirely by assuming that the tangential velocity of the interface was distinct from the underlying potential 
flow. This tangential velocity could be calculated by integrating the second-order differential equation 
that resulted from the tangential stress boundary condition. This second-order equation was solved by 
representing the tangential velocity of the interface as a truncated sine series and solving for the coefficients 
in a least-squares sense. Once the surface tangential velocity was known, the surfactant transport equation 
could be used to update the surfactant concentration simultaneously with the remaining surface variables 
and the drop shape. The resulting calculations for effects of surface viscosity were again found to agree with 
small-amplitude theory. 
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Appendix A 

Surface theorems 


Here the general forms for the Surface Reynolds Transport Theorem and Surface Divergence Theorem are 
given. The reader is referred to [4, 16, 31] for their derivation. 

The Surface Reynolds Transport Theorem for a material surface S m is: 


d 

dt JSi 


L* dS =L 


~dt 


+ V s • (v#) dS . 


(A.l) 


where is a scalar, vector, or tensor field. 

Referring to Figure A.l, the surface Divergence Theorem for a material surface S rn and material line 
element C m is: 

r r 

(A. 2) 


/ V,-(I f D*)dS= / Tja&dl, 
Js m Jc m 


where x ”, or “(blank)”. 



Figure A.l: A material surface S m bounded by the material line C m , where f) and f are perpendicular to 
the unit normal n. 
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Appendix B 

Constants in matched asymptotics 


Here the constants in the time-periodic uniformly valid velocity profiles and surfactant concentrations are 
given explicitly for cases 1-3. 


B.l Case 1: “negligible” surfactant effects 

The constants Cq 0 and Cq 0 are 

C — - 11 vgg 

60 “ (t + 1) (y/iip + Vw) 

- (2£+l) y/w 

60 t (x/W + y/£p) ' 


The constants Cex and Cox are given by the solution of the simultaneous equations: 


—ICqi — (£ + \)C$\ — — 


-fCoo + (£+ 1 fC do 
a 


£Cq i + (l + l)Coi — — l —£[{(' ~h 2) + Cq 0 \ — (£ -h !)[(£ — 1) Cg 


j} 


B.2 Case 2: “small” surfactant effects 

The constants C$ i and Ce 2 are 

-1 y/2 


Cox = 
C 02 = 


eo (1 + 0 

1 V2 


V^eo (1 + 0 

The constants G, C g \ . and C g 2 are 


2(1 - 1) + %g + ^(e - i)(f + 2) + ^e.(e + 1) 

or nr 

-Cex \2-%%±-^(e-m + 2)-^e(£+ 1) 

Q 2 Cox (X a 


G = £(£ — l)Pe, 
(f + 1) 


t(t + 1) — tfhoPc* 


Cg 1 = 

Cg2 = 


(f-1) 

(<+l) 

(f-1) 


Cex 

C 9 2- 


[p(e + i)i + n> 0 Pe 2 s \ 

V 2 


(i + O 


' + 1) 


Cex 

\fiho 


(B.l) 

(B.2) 

(B.3) 

(B.4) 


(B.5) 

(B.6) 


(B.7) 

(B.8) 

(B.9) 
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B.3 Case 3: “medium” surfactant effects 


The constant Cq o is 


where 


The constant G is 


The constants Cq 
They are not needed 


_ ni + n 2 e l7r/2 
00 “ d x + d 2 e^/ 4 + d 3 e™/ 2 5 


(B.10) 


n i = - 


^2 

c/i 

d‘2 


e*Pe„l 2 (l - l)(l + 1) 
*(* + l) 2 + n 2 0 Pe 2 
e:n«Pe^(£-l) 


-i)(<+ 2) 


^ + l) 2 + f^ 0 Pe; 

e*Pe s f 2 (£ + l) 2 


f 2 (£+l) 2 +n 2 0 Pe 

QV^fO 

6sf^oPe 2 ^(^ + 1) 
~l?\l + l) 2 + fl 2 0 Pe; 


? + K ;^ + i) + /i;(^- 1)^ + 2) 


(B.ll) 

(B.12) 

(B.13) 

(B.14) 

(B.15) 


G - l{l- l)Pe, 


t{JL + 1 ) - iShoPes ' 

e*(i + i ) 2 + u 2 0 p e 2 . ' 


(B.16) 


and C g i take similar but more complicated forms and are not written explicitly here, 
in the leading-order analysis of the energy equation for Case 3. 
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Appendix C 

Local coordinate system 


x(s,<£>) = e r (p)r(s) + e 2 z(s) (C.l) 

s (s,p) = e r (y>)r'(s) + e 2 z'(s) (C.2) 

n(s,<^) = — e r (c^)z / (5) + e 2 r'(s) (C.3) 

In local coordinates (n,s,c^), the unit tensor, the gradient, operator and their projections onto the interface 
are given by 

I = nn + ss + , (C.4) 

I 5 = I - nn = ss + , (C.5) 

_ A , d A d . , d (n 

V — n h n — ~h s h s ~ h ^ — , ( L- . 0 ) 

dn ds dp 

d d 

V s = I s • V = sh s — + e* /i^— , (C.7) 

where h n ,h s , and h 9 are metric coefficients that are functions of n,s, and p in general. For points at a 
constant meridian angle ( p = constant) and on the interface (n = 0), the metrical coefficients and their 
derivatives are given by 


h n = 1 


h s — 1 


= - 
r 


i- 

dn _ 
rj — S C S 
as 


= — n C, 


^ = ^ =nZ ~ Sr ’ 

where C s = z'r" - r'z" is the curvature in the s direction, and primes denote differentiation with respect to 
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Appendix D 

Explicit kernels in numerics 


Here the forms of the kernels and their asymptotic behavior as s -> Si are given explicitly in terms of the 
axisymmetric arclength parameterization (at constant <p) x(s) = e r r(s) + e z z(s), where ds 2 = dx • dx. The 
kernels in the regularized double-layer potential boundary integral representation (3.8) and relation (3.9) 
take the parameterized forms 

h(s) * VG(sj,s)r(s) - 

r(s){[z(sj) - z(s)]r ; (s) + [r(sj) + r(s)]z'(s)} f E{m) 

7r{[r(si) + r(s)] 2 + [z(si) - z(s)] 2 } 3 / 2 [(1 - m) 

+ - K(m) } , (D.l) 

27r{[r(sj) + r(s)] 2 + [*(*) - z(s)] 2 }!/2 [(1 - m) 'J 


where 


e ¥ -[n(s)x VG(Si,s)r(s)| = 

r(s){[z(sj) - z(s)]z'(s) - \rjsj) + r(s)]r'(s)} Ejm ) 

7r{[r(sj) + r(s)] 2 + [z(si) - z(s)] 2 } 3 / 2 [(1 - m) 

{[ 2 {s;) - z(s)]z'(s) - r(s)r'(s)} [ Ejm) 


27rr(si){[r(s t ) + r(s)] 2 + [z(s;) -z(s)] 2 } 1 / 2 [(1 -m) 


— K(m) 


, v 4r(sj)r(s) , 

m m(si,s) [r(s . )+r(s)] 2 + Ws .)_ 2(s )]2 ‘ 

and primes (/) denote differentiation with respect to arclength s . K(m) and E(m) are, respectively, the 
complete elliptic integrals of the first and second kind. Accurate asymptotic expressions for K(m) and E(m) 
may be found in [1]. 

As s — > Si the above kernels have the asymptotic behaviors 
h(s)VG(s u s)r(s) ~ 

T(-C y ( B< )logf , (S ~ Si) . 1 +C,(a i )-C v (s i )} + 0(a-s i ), P-4) 

4tt l. L4\/2r(s i )J J 


b v ■ [n(s) x VG(sj, s) r(s)] ~ 

— J 1 . I r '( s ») lo [ ( g " s ») 

27T \ (s - Si) 2 r(si) _4y/2r(si) 


+ G(s — si) , 


where C s (s) = z'(s)r"(s) - r'(s)z"(s) and C v {s) = -z'(s)/r(s). 
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